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1.  INTRODUCTION 

1.1  General 

Various  engineering  problems  can  be  reduced  to  the  solution  of  matrix 
eigenvalue  problems.  Typical  examples  in  the  field  of  structural  engineering 
are  the  problem  of  determination  of  natural  frequencies  and  the  corresponding 
normal  modes  in  a  dynamic  analysis  and  the  problem  of  finding  buckling  loads 
in  a  stability  analysis  of  structures.  Since  the  advent  of  the  digital  com¬ 
puter,  the  complexity  of  structures  which  can  be  treated  and  the  order  of 
the  corresponding  eigenvalue  problems  have  been  greatly  increased.  Hence, 
the  development  of  solution  techniques  for  such  problems  has  attracted  much 
attention. 

For  the  dynamic  analysis  of  a  linear  discrete  structural  system  by 
superposition  of  modes,  we  must  first  solve  the  problem  of  free  vibration 
of  the  system.  The  free  vibration  analysis  of  the  linear  system  without 
damping  reduces  to  the  solution  of  the  linear  eigenvalue  problem 

Ax  =  X  3x  (1.1) 

in  which  A  and  B  are  stiffness  and  mass  matrices  of  order  n,  the  number  of 
degrees  of  freedom  of  the  structural  system.  A  column  vector  x  is  an 
eigenvector  (or  normal  mode),  and  the  scalar  x  the  corresponding  eigenvalue 
(or  the  square  of  a  natural  frequency). 

The  matrices  A  and  B  are  real  and  symmetric,  and  are  usually  banded  and 
sparse.  If  a  consistent  mass  matrix  is  used,  the  matrices  A  and  B  have  the 
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same  bandwidth  [4,5].  If  a  lumped  mass  model  of  the  system  is  used,  B  will 
be  diagonal.  The  matrix  B  is  positive  definite,  but  the  matrix  A  may  be 
semidefinite.  There  are  n  sets  of  solutions  of  Eq.  (1.1),  that  is,  n  eigen¬ 
values  and  their  corresponding  eigenvectors. 

Frequently,  in  practical  eigenvalue  problems,  the  order  of  A  and  B  is 
so  high  that  it  is  impractical  or  very  expensive  to  obtain  the  complete 
eigensolution.  On  the  other  hand,  to  carry  out  a  reasonably  accurate 
dynamic  analysis  of  the  structure,  it  is  possible  to  consider  only  a  partial 
eigensolution.  The  partial  solution  of  interest  may  consist  of  only  few 
lowest  eigenvalues  and  their  eigenvectors,  or  eigenvalues  in  the  vicinity 
of  a  given  frequency  and  the  corresponding  eigenvectors.  The  method 
described  in  this  study  is  aimed  at  effective  solution  of  this  type  of 
problem  rather  than  at  a  complete  eigensolution. 

1.2  Object  and  Scope 

The  object  of  this  study  is  to  present  an  iterative  method  which  is 
efficient  and  numerically  stable  for  the  accurate  computation  of  limited 
number  of  eigenvalues  and  the  corresponding  eigenvectors  of  linear  eigenvalue 
problems  of  large  order. 

The  method  developed  remedies  the  major  drawbacks  of  the  inverse  iter¬ 
ation  method  with  spectral  shifting  [13]:  numerical  instability  due  to 
shifting  and  slow  convergence  when  eigenvalues  are  equal  or  close  in  magni¬ 
tude.  The  proposed  method  converges  rapidly  and  is  numerically  stable  for  any 
number  of  multiple  or  close  eigenvalues  and  the  corresponding  eigenvectors. 
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The  procedure  for  distinct  eigenvalues  is  treated  in  Chapter  2,  and  a 
modified  procedure  for  multiple  or  close  eigenvalues  in  Chapter  3.  Selection 
of  initial  approximate  eigenvalues  and  eigenvectors  by  the  subspace  iteration 
method  is  described  in  Chapter  4.  To  show  the  efficiency  of  the  proposed 
method,  three  sample  problems  are  solved:  vibration  of  a  plane  frame,  of 
a  plate  in  bending,  and  of  an  arch.  Comparisons  are  made  in  Chapter  5  with 
a  method  which  is  generally  regarded  as  very  efficient,  the  subspace 
iteration  method. 

1.3  Review  of  Solution  Methods 

Numerous  techniques  for  the  solution  of  eigenvalue  problems  have  been 
developed.  These  techniques  can  be  divided  into  two  classes  -  techniques 
for  approximate  solution  and  techniques  for  "exact"  solution. 

The  approximate  solution  techniques  include  well-known  static  conden¬ 
sation  [2,3,24,25,27,42],  dynamic  condensation  [34],  Rayleigh-Ritz  analysis 
[9,13,31,48],  component  mode  analysis  and  related  methods  summarized  by 
Uhrig  [50].  These  methods  are  essentially  techniques  for  reducing  the  size 
of  a  system  of  equations.  The  reduction  of  a  system  of  equations  eventually 
leads  to  a  loss  in  accuracy  of  a  solution.  However,  the  advantage  of 
lessened  computational  effort  for  a  solution  sometimes  may  compensate  for 
the  loss  in  accuracy.  Moreover,  an  approximate  solution  found  by  these  methods 
may  serve  as  the  starting  solution  for  the  exact  methods,  which  will  be 
discussed  next. 

The  exact  methods  are  designed  for  the  accurate  computation  of  some 
or  all  the  eigenvalues  and  corresponding  eigenvectors.  These  methods  consist 
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of  vector  iteration  methods,  transformation  methods,  the  method  based  on 
the  Sturm-sequence  property,  polynomial  iteration  method,  and  minimization 
methods.  These  methods  are  well  described  in  Ref.  51.  The  methods  differ 
in  the  choice  of  which  mathematical  properties  of  an  eigenvalue  problem  are 
used.  The  vector  iteration  methods  such  as  the  classical  vector  iteration 
(power  method)  and  simultaneous  vector  iteration  deal  with  the  form  of 
equations  Ax  =  x  Bx.  The  transformation  methods  (LR,  QR,  Jacobi,  Givens, 
and  Householder  methods)  are  based  on  the  mathematical  property  that  the 
eigenvalues  of  a  system  are  invariant  under  similarity  transformations.  In 
the  polynomial  iteration  method,  the  roots  of  det  (A  -  xB  )  =  0  are  found, 
and  minimization  methods  are  based  on  the  stationary  property  of  the 
Rayleigh  quotient  [43]. 

In  vector  iteration  methods  and  minimization  methods,  both  the  eigen¬ 
values  and  corresponding  eigenvectors  are  found  simultaneously,  but  in 
other  exact  methods,  only  eigenvalues  are  computed  or  the  computed  eigen¬ 
vectors  are,  in  general,  not  suitable  for  use  in  the  final  solutions.  In 
such  methods,  another  method  such  as  the  vector  iteration  method  with  a 
shift  may  be  used  for  finding  the  eigenvector  corresponding  to  a  computed 
eigen val ue. 

For  a  limited  number  of  eigenvalues  and  corresponding  eigenvectors  of 
an  eigenvalue  problem  of  large  order  which  we  are  concerned  with  in  this 
study,  the  above  methods  have  been  modified  or  combined  to  take  advantage 
of  the  useful  characteristics  of  several  of  the  methods.  First,  the 
determinant  search  method  [7,9,22,23]  combines  the  methods  based  on  the  Sturm- 
sequence  property,  polynomial  iteration,  and  inverse  iteration.  In  this 
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method,  eigenvalues  in  a  specified  range  are  approximately  isolated  by  using 
the  bisection  method  and  the  Sturm-sequence  property  and  then  located 
accurately  by  the  polynomial  iteration  method.  The  corresponding  eigen¬ 
vectors  are  computed  by  inverse  iteration  with  a  shift.  By  this  method, 
eigenvalues  in  any  range  and  corresponding  eigenvectors  can  be  found. 

However,  it  has  the  disadvantage  that  the  matrix  is  factorized  in  each  iter¬ 
ation  to  locate  the  eigenvalues  of  interest. 

Another  method  for  the  solution  of  large  eigenvalue  problems  is  the 
so-called  subspace  iteration  method  [6,15,32,39,47],  which  is  a  combination 
of  the  simultaneous  iteration  method  and  a  Rayleigh-Ritz  analysis,  in  this 
method,  several  independent  vectors  are  improved  by  vector  inverse  iteration, 
and  the  best  approximation  to  the  eigenvectors  are  found  in  the  subspace  of 
the  iteration  vectors  by  a  Rayleigh-Ritz  analysis.  In  this  method,  eigen¬ 
values  at  the  end  of  the  spectrum  and  the  corresponding  eigenvectors  converge 
very  rapidly.  This  method  will  be  discussed  further  in  Chapter  4. 

The  inverse  iteration  method  with  a  shift  is  known  to  be  extremely 
efficient  for  improving  approximate  eigenvalues  and  eigenvectors.  However, 
as  mentioned  in  the  previous  section,  when  the  shift  is  very  close  to  a  true 
eigenvalue,  the  method  exhibits  numerical  instability,  yielding  unreliable 
answers  [13].  In  addition,  when  the  eigenvalues  of  interest  are  close  to¬ 
gether,  their  convergence  is  very  slow.  Robinson  and  Harris  [44]  developed 
an  efficient  method  to  overcome  the  above  difficulty  for  distinct  eigenvalues 
by  augmenting  the  equations  used  in  the  inverse  iteration  method  by  a  side 
equation.  While  this  method  extracts  eigenvalues  and  eigenvectors  simul¬ 
taneously  with  a  very  high  convergence  rate,  it  has  the  disadvantage  that  the 
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algorithm  is  inefficient  for  problems  with  multiple  or  close  eigenvalues. 
This  method  and  some  improvements  on  it  will  be  discussed  further  in  the 
next  chapter. 


1.4  Notation 

All  symbols  are  defined  in  the  text  when  they  first  appear. 

With  regard  to  matrices,  vectors,  elements  of  matrices  or  vectors,  and 
iteration  steps,  the  following  conventions  are  generally  used: 

(1)  Matrices  are  denoted  by  uppercase  letters,  as  A,  B  and  X. 

(2)  A  column  vector  is  denoted  by  a  lowercase  letter  with  a 
superior  bar  and  a  subscript,  as  a.,  5.  and  x.. 

J  J  J 

(3)  Elements  of  a  matrix  or  vector  are  denoted  by  a  lowercase 
letter  with  a  double  subscripts,  as  a,,,  b. .  and  x. .. 

I J  I J  I J 

-(k\ 

(4)  Iteration  steps  are  denoted  by  a  superscript,  as  X'  *  xj 
and  x{« 

fkl  (k) 

(5)  Increments  are  denoted  by  the  symbol  A,  as  ax^  '  and  Axi.. 

Some  symbols  are  assigned  more  than  one  meaning.  However,  in  the  context 
of  their  use  there  are  no  ambiguities. 


A,  Sj,  a(j 
A*(k) 


B,  bj.  blJ 


,*(k) 


XL 

=  stiffness  matrix,  j u  column  vector  of  A,  element 
of  A 

=  projection  of  A  onto  the  subspace  spanned  by  vectors 
in  Y^,  A*<k)  =  Y*k)T  A 

=  radius  of  circular  arch 

!■  L. 

=  mass  matrix,  jin  column  vector  of  B,  element  of  B 

=  projection  of  B  onto  the  subspace  spanned  by  vectors 
in  Y^k\  B*^  =  Y^T  B  Y^k) 
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expansion  matrix  of  X^,  j**1  column  vector  of  C^, 

element  of  C^,  X^  =  XC^ 

diagonal  matrix,  see  Section  2.2 

plate  bending  stiffness,  Dfi  =  EH^/12(l-u^) 

matrix  for  finding  close  or  multiple  eigenvalues  and 
eigenvectors,  jth  column  vector  of  D,  see  Eq.  (3.24) 
iteration  matrix  for  D  after  k  iterations,  jth  column 
vector  of  D^,  see  Eq.  (3.23) 

Young's  modulus 

diagonal  matrix,  jth  column  vector  of  E,  element  of 
E,  see  Eq.  (A. 7) 

★ 

diagonal  matrix,  jth  column  vector  of  E  ,  element  of 


thickness  of  plate 

number  indicating  rate  of  convergence  of  eigenvector, 
see  Eq.  (2.13) 

moment  of  inertia  of  cross-section 
identity  matrix  of  order  s 

indices  of  matrix  elements 

superscript  indicating  number  of  iterations 

lower  triangular  matrix 

Lagrangian,  see  Eq.  (3.6) 

average  half  bandwidth  of  A,  of  B 

total  number  of  operations  required  for  finding 
eigenpairs  by  the  proposed  method,  by  the  Robinson- 
Harris  method,  by  the  subspace-iteration  method 
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n 

P 

q 


s 


X*  xj*  xij 


x(k)  =(k)  x(k) 

X  ,  Xj  ,  x^ 


vOO  y(k)  v(k) 
Y  ’  yj  *  yij 


Z,  Z 
Y(k) 


(k) 


6ij 

e(k) 

6j 

* 

x 

A,  A, 


order  of  A  and  B 

number  of  eigenpairs  sought 

number  of  iteration  vectors  by  subspace  iteration 

method,  q  =  max(2p,  p+8) 

residual  vector  of  approximation  to  jth  eigenpair 
after  k  iterations 

number  of  close  and/or  multiple  eigenpairs  sought 
number  of  iterations  needed  to  find  eigenpairs  by 
proposed  method,  by  Robinson-Harris  method,  by 
subspace  iteration  method 

matrix  of  eigenvectors  (modal  matrix),  jth  eigen¬ 
vector,  element  of  X 

approximation  to  X  after  k  iterations,  jth  column 
vector  of  X^,  element  of 

(kl 

matrix  of  iteration  vectors  improved  from  Xv  '  by 
simultaneous  iteration  method,  jth  column  vector  of 
element  of 

rotation  matrix,  approximation  to  Z  after  k  iterations 
error  in  x(^  or 

increment  operator 
Kronecker  delta 

error  in  x^)  or  yi^ 

multiple  eigenvalue 

diagonal  matrix  of  eigenvalues,  jth  eigenvalue, 

A  =  diag(X-| ,  Xg,  •  •  •  , 
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A(k)  x(k) 

’  j 


V 

p1j* 


p 


0) 


=  approximation  to  A,  to  a.,  after  k  iterations 

J 

*  shift  applied  in  vector  iteration  method 

*  element  of  D,  of 

=  mass  density 

2 

=  natural  circular  frequency,  x  =  w 
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2.  DISTINCT  ROOTS 


2.1  General 

In  this  chapter,  a  method  for  finding  a  simple  eigenvalue  and  the 
corresponding  eigenvector  will  be  presented.  The  method  developed  by 
Robinson  and  Harris  [44]  is  modified  here  to  save  overall  computational 
effort  for  finding  an  eigensolution.  The  Robinson-Harris  method  is  an 
application  of  the  Newton-Raphson  technique  for  improving  the  accuracy  of 
an  approximate  eigenvalue  and  the  corresponding  approximate  eigenvector. 

In  the  proposed  method,  a  modified  form  of  the  Newton-Raphson  technique  Is 
applied  instead  of  the  standard  one  used  in  the  Robinson-Harris  method. 

In  Section  2.2,  the  Robinson-Harris  method  will  be  discussed  first; 
then  the  proposed  method  will  be  presented.  The  convergence  rate  of  the 
proposed  method  and  the  number  of  operations  per  iteration  will  be  given  in 
Section  2.3.  The  estimation  of  error  in  an  approximate  solution  is  found 
in  Section  2.4.  A  technique  for  the  examination  of  the  converged  solution 
to  determine  whether  the  eigenvalues  and  corresponding  eigenvectors  of 
interest  have  been  missed  and  a  method  for  finding  a  missed  solution  will 
be  presented  in  Section  2.5. 

2.2  The  Iterative  Scheme 

Let  us  consider  the  following  linear  eigenvalue  problem 

AXj  =  Aj  BXj  (j  *  1,  2,  ...»  n)  (2.1) 


I 
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where  A  and  B  are  assumed  to  be  given  symmetric  matrices  of  order  n  and  B 

Is  taken  to  be  positive  definite.  The  x.  and  x.  are  the  jth  eigenvalue  and 

J  J 

the  corresponding  eigenvector. 

Let  us  assume  that  an  initial  approximate  solution  of  Eq.  (2.1),  x.^ 
and  is  available.  Denote  an  approximate  eigenvalue  and  the  corre¬ 

sponding  eigenvector  after  k  iterations  by  x.^  and  x.^  (k  =  0,  l,  .  .  .) 

J  J 

Then,  we  have 


Ax.(k)  -  x<k>  Bx.^  =  r.^ 

J  J  J  J 


(2.2) 


-  (L) 

where  r.'  '  is  a  residual  vector. 

v 

The  object  is  to  remove  the  residual  vector  in  Eq.  (2.2).  The  Newton 
Raphson  technique  is  applied  for  this  purpose.  Let  the  (k  +  1 ) th  approxi¬ 
mation  be  defined  by 


(k+1)  _  .  (k) 
j  "  Aj 


(k) 


-  (k+1)  -  (k)  -  (k) 

xj  '  *j  +  AXj 


(2.3) 


(k)  -  (k)  (k) 

where  ax.'  and  ax-'  '  are  small  unknown  incremental  changes  of  x.'  '  and 

w  J  J 

S<k>.  Substituting  X-^k+^  and  x-^k+1^  of  Eq.  (2.3)  for  and  x.  in 

J  J  J  J  J 

Eq.  (2.1)  and  discarding  a  nonlinear  term  ax  BAx-^k^  as  very  small 

J 

compared  with  the  other,  linear,  terms,  we  get 


(A  -  X.(k)B)  Ax.(k)  -  AX.(k)  Bx.(k)  =  -  (2.4) 

J  J  J  j  J 


where  r/k^  is  the  residual  vector  defined  in  Eq.  (2.2). 
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Note  that  in  Eq.  (2.4),  there  are  n+1  scalar  unknowns  (aa.^  and  n 

-  (k)  J 

components  of  ax.  '),  but  only  n  equations.  Hence,  it  is  required  for  the 

J 

solution  of  Eq.  (2.4)  that  either  the  number  of  unknowns  be  reduced  or  one 

equation  added.  Derwidue  [16]  and  Rail  [41]  reduced  the  number  of  unknowns 

by  setting  the  nth  component  of  the  vector  ax.^  or  x.^k+1^  at  a  preassigned 

J  J 

value  -  zero  or  one.  In  these  methods,  it  may  happen  that  an  unfortunate 
choice  of  one  component  results  in  failure  of  the  procedure. 

Instead  of  reducing  the  number  of  unknowns,  Robinson  and  Harris  [44] 
added  an  extra  equation  (side  condition)  to  the  system  of  Eq.  (2.4),  to 
arrive  at  a  set  of  n+1  equations  in  n+1  unknowns.  This  side  condition  is 

-x<k>T  BAx.(k)  =  0  (2.5) 

J  J 

_  / 

Equation  (2.5)  means  that  the  incremental  value  ax-  '  is  orthogonal  to 

J 

-  Ik) 

the  current  approximate  eigenvector  x.'  '  with  respect  to  the  matrix  B. 

J 

-  (k) 

The  side  condition  prevents  unlimited  change  in  the  x.v  .  The  resulting 

J 

set  of  simultaneous  linear  equations  may  be  written  in  matrix  form  as 


A  -  a.^B 

J 

1 

-  B5j(k) 

1  X 

<1 

1 _ 

II 

: 

o 

1 _ 

1 

> 

C_i. 

7T 

0 

■  i 

_  /i.\ 

where  the  residual  vector  r.K  1  is  given  in  Eq.  (2.2).  The  coefficient 

J 

matrix  for  the  incremental  values  is  of  order  n+1  and  symmetric.  Moreover, 

it  is  nonsingular  if  Aj  is  not  multiple  [44].  Equation  (2.6)  may  be  solved 
(k)  -  (k) 

for  AAj  '  and  AXj  '  by  Gauss  elimination,  or  by  any  other  suitable 


i 


1 


m 
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technique.  Note  that  the  submatrix  in  the  coefficient  matrix  (A  -  a.^B) 

(k) 

is  almost  singular  when  A.v  '  is  close  to  a..  However,  this  does  not  cause 

J  J 

any  difficulty  in  solving  Eq.  (2.6),  since  in  the  elimination  process  only 
the  last  pivot  element,  in  general,  becomes  very  small.  Thus,  the  inter¬ 
change  of  columns  and  rows  does  not  increase  significantly  the  column  height 

of  the  factorized  matrix.  The  improved  values,  A.^k+^  and  x.^+^,  are 

J  3 

computed  from  Eq.  (2.3).  The  procedure  employing  Eqs.  (2.3)  and  (2.6)  is 

repeated  until  the  errors  in  the  A.^)  and  x  ^  are  within  allowable  toler- 

J  J 

ances.  The  method  of  estimating  these  errors  will  be  discussed  in  Sec¬ 
tion  2.4. 

The  convergence  of  the  above  process  for  an  eigenvalue  and  the  corre¬ 
sponding  eigenvector  has  been  shown  to  be  better  than  second  order;  the 
order  has  been  found  to  be  2.41  [44].  However,  the  algorithm  using  Eq.  (2.6) 
requires  a  new  triangularization  in  each  iteration,  since  the  values  of  the 
elements  of  the  coefficient  matrix  are  changed  in  each  iteration  as  a  result 
of  changing  from  A.  ^  to  A.^k+1^.  The  number  of  operations  (multiplications 

J  J 

and  divisions)  required  in  such  a  triangularization  is  very  large. 

To  avoid  the  complete  elimination  procedure  in  each  iteration,  the 
following  equations  instead  of  Eq.  (2.6)  are  used  in  the  proposed  method. 


“ 

• 

- 

A  -  A.{0)B 

J 

4;<k> 

sj 

-  5.(k)  B 

m  J 

0 

- 

*< 

< 

_ 1 

1 

o 

_ > 

(2.7) 


GO 
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where  the  residual  vector  rS^  is  defined  in  Eq.  (2.2).  Equation  (2.7) 
was  obtained  by  introducing  Eq.  (2.3)  into  Eq.  (2.1)  and  discarding  a  small 
linear  term  (x.^k+1^  -  a.^)  Bax.^.  Note  that  Eq.  (2.7)  differs  from 

J  v  J 

Eq.  (2.6)  in  such  a  way  that  the  coefficient  matrix  in  Eq.  (2.6)  has  the 
submatrix  (A  -  a.^B),  while  the  coefficient  matrix  in  Eq.  (2.7)  has 

J 

(A  -  A.^B).  The  coefficient  matrix  in  Eq.  (2.7)  is  also  symmetric,  and 

J 

nonsingular  if  A.  is  not  multiple.  The  nonsingularity  of  the  coefficient 

J 

matrix  will  be  proved,  in  passing,  in  Appendix  A. 

From  the  form  of  the  coefficient  matrix,  it  can  be  seen  that  once  the 
matrix  is  decomposed  into  the  form  LOlJ,  where  L  is  lower  triangular  and  D 
is  diagonal,  only  a  small  number  of  additional  operations  is  required 
for  the  solution  of  Eq.  (2.7)  in  the  succeeding  iterations,  since  only  the 
vector  Bx.^k)  in  the  matrix  is  changed  in  each  iteration.  The  proposed 

J 

method  therefore  considerably  reduces  the  number  of  operations  required  in 

each  iteration.  On  the  other  hand,  the  method  lowers  the  convergence  rate 

because  of  the  neglect  of  the  small  linear  term  (x-^k+1^  -  a.^)  (Bax-^), 

J  J  J 

which  in  turn  increases  the  number  of  iterations  for  a  solution.  However, 
the  overall  computational  effort  for  a  solution  does  decrease.  It  will  be 
seen  in  Chapter  5  that  the  proposed  method  is  actually  more  efficient  than 
the  Robinson-Harris  method. 

2.3  Convergence  Rate  and  Operation  Count 

The  efficiency  of  a  numerical  method  such  as  the  one  proposed  here  can 
be  estimated  given  the  convergence  rate  and  the  number  of  operations  per 
iteration  required  in  the  process.  The  convergence  analysis,  which  is  given 


** 
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in  Appendix  B,  will  be  summarized  as  follows.  Let  an  approximate  eigenvector 
-  (k) 

x.'  be  expanded  in  terms  of  the  true  eigenvectors  x-,  i.e., 

J  • 


-  (k)  =  -  c  (k)  - 
xj  -  ci  j  xi 

i  =  l 


(2.8) 


(kl  -  fkl 

where  c.jK  '  is  a  coefficient  of  the  vector  x^.  If  Yj  is  the  error  in 

and  0.^  the  error  in  x.^,  they  may  be  defined  as 

J  J  v 


(k)  _ 


A.  -  A. 

_ 1 


(k) 


Xj 


(2.9) 


e. 

J 


(k)  _ 


i=l 

l&L 


i=l 


1/2 


(2.10) 


where  0 


<*> , 


s  a  measure  of  the  angle  between  the  vectors  c 


(k) 


-  fkl 1 

where  c.v  7  * 

J 


(CU 


(k) 


*2j 


(k) 


and  Cj*  and 

cn j  and  cj  =  cjj{k).0,...»0). 


The  geometric  interpretation  of  e.^  is  illustrated  in  Fig.  1. 

J 

With  the  above  definitions,  the  errors  in  A.^+1)  and  x.^k+1^  may  be 

J  J 

written  as  (see  Appendix  B) 


T,(W>  *  T  <« 


(2.11) 


1 
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e.(k+1)  =  he.(k)  (2.12) 

J  J 

where 

A.  -  X.<°> 

h  =  -J - J.  T  <  1  (m  =  1 ,  2,  .  .  .  n)  (2.13) 

A  -  A  S°>  ~ 

m  j 

Equations  (2.11)  and  (2.12)  show  that  the  convergence  character  of  both 
eigenvalues  and  eigenvectors  is  linear.  However,  the  eigenvalues  converge 
much  more  rapidly  than  the  eigenvectors.  Note  also  that  the  closer  a^.  is 
to  another  eigenvalue,  the  larger  a  is,  yielding  slow  convergence.  Hence, 
the  method  is  not  suitable  for  finding  close  eigenvalues  and  the  corre¬ 
sponding  eigenvectors. 

Another  important  consideration  which  should  be  taken  into  account  in 
estimating  the  efficiency  of  numerical  methods  is  the  number  of  operations 
per  iteration.  One  operation  is  defined  as  one  multiplication  or  division, 
which  almost  always  is  followed  by  an  addition  or  a  subtraction.  For  the 
expression  of  this  number,  let  m3  and  m.  be  the  half  band-widths  of  the 
matrices  A  and  B,  and  let  n  be  the  order  of  A  and  B.  Let  Tp  be  the  number 
of  iterations  needed  to  find  p  eigenpairs  by  the  proposed  method  and  Tr  by 
the  Robinson-Harris  method.  Then,  the  number  of  operations  for  p  eigenpairs, 
Np,  required  by  the  proposed  method  is 


N  =  -^pn  (m2  +  3m  +  2m.  +  2)  +  Tn  (5ma  +  2m.  +  6) 
p  2  a  a  b'p  a  b 


(2.14) 


*  » 


I 

I 

1 1 


and  by  the  Robinson-Harris  method,  Nr,  is 


Nr  =  iTrn  ^ma  +  13ma  +  6mb  +  12  ^ 


(2.15) 


It  can  be  seen  that  the  number  of  operations  per  iteration  required  by  the 
proposed  method  is  much  smaller  than  for  the  Robinson-Harris  method.  The 
development  of  the  above  expressions  is  given  in  Table  1. 


2.4  Errors  in  Approximate  Eigensolutions 


An  important  feature  of  an  iterative  method  such  as  the  proposed  method 

is  some  means  of  estimating  the  error  in  a  computed  solution.  This  permits 

one  to  terminate  the  iteration  process  at  the  point  where  a  sufficiently 

accurate  result  has  been  obtained.  It  is  important  to  have  estimates  in 

terms  of  numbers  available  in  the  calculation,  since  it  is  impossible  to 

compare  with  the  exact  values. 

(kl  ( kl 

The  error  in  x.v  ,  y- '  ,  can  be  estimated  as  follows:  from 

J  J 

Eqs.  (2.9)  and  (2.11) 


(2.16) 


Substituting  Eq.  (2.16)  for  A.  in  Eq.  (2.9)  gives 

J 


M  . 


,  (k+l)  fl  ±  ,2  T(k)  ,  (k+l) 

J  J  J 


(2.17) 
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Since  0  <  h  «  1  and  0  <  y.: 

”  J 


«  1 ,  from  Eq.  (2.17) 


The  error  in  x/k\  e/k\  can  be  approximated  by  [e.^  -  e.^k+1^] 
J  J  j  j 

since  e/k+^  «e.^.  Furthermore,  from  Fig.  1, 

J  J 


(2 


(2 


filial 
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The  number  of  operations  for  the  estimation  of  e.^  is  only  about 

J 

n  (2m^  +  3),  which  is  small  compared  with  the  number  of  operations  per  iter¬ 
ation  (see  Section  2.3). 

2.5  Treatment  of  Missed  Eigensolutions 

Some  of  the  eigenvalues  and  corresponding  eigenvectors  of  interest  may 

be  missed  when  the  initial  approximations  are  not  suitable.  In  order  to 

check  whether  this  occurs,  the  Sturm-sequence  property  [9,31,39,48,51]  may 

be  applied.  The  Sturm-sequence  property  is  expressed  as  follows:  if  for 

an  approximate  eigenvalue  A.^,  (A  -  A.^B)  is  decomposed  into  LDLT, 

J  J 

where  L  is  a  lower  triangular  matrix  and  D  a  diagonal  one,  then  the  number 
of  negative  elements  in  D  equals  the  number  of  eigenvalues  smaller  than 
A.^.  A  computed  eigenvalue  can  be  checked  using  the  above  property  with 

J 

negligible  extra  computation,  since  the  decomposition  of  the  matrix 
(A  -  Aj^B)  has  already  been  carried  out  during  the  procedure  for  the 
solution  of  Eq.  (2.7). 

If  some  of  the  eigenvalues  of  interest  are  detected  to  be  missing, 
finding  them  consists  of  three  steps:  finding  approximations  to  the  mis;ed 
eigenvalues,  finding  approximate  eigenvectors  corresponding  to  the  missed 
eigenvalues,  and  improving  the  approximate  eigensolutions. 

The  approximate  eigenvalues  can  be  found  by  the  repeated  applications 
of  the  Sturm-sequence  calculation  mentioned  above  and  the  method  of  bisection 
[9,31,38,51],  or  by  the  polynomial  iteration  method  [7,8,9,38,51],  in  which 
the  zeros  of  the  characteristics  polynomial  p(A)  =  det(A  -  AB)  are  found 
using  variants  of  Newton's  method. 


In  the  second  step,  the  approximation  to  the  eigenvectors  corresponding 
to  the  missed  eigenvalues  is  found.  Frequently,  finding  the  eigenvectors 
corresponding  to  the  missed  eigenvalues  is  much  more  difficult  than  finding 
the  missed  eigenvalues.  However,  subspace  iterations  with  a  shift  [6,32], 
which  will  be  discussed  in  Chapter  4,  or  dynamic  condensation  [34,42,50] 
may  be  used  for  this  purpose. 

Finally,  the  approximate  eigenvalues  and  corresponding  approximate 
eigenvectors  can  be  improved  by  the  method  of  Section  2.2  if  the  eigen¬ 
values  are  not  multiple  or  close,  or  if  they  are,  by  the  method  of 
Chapter  3. 
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3.  CLOSE  OR  MULTIPLE  ROOTS 


3. 1  General 

As  mentioned  earlier,  the  method  presented  in  Chapter  2  fails  or 
exhibits  slow  convergence  if  it  is  applied  to  the  solution  for  multiple  or 
close  eigenvalues  and  for  their  corresponding  eigenvectors.  The  failure  or 
slow  convergence  of  the  method  is  caused  by  impending  singularity  of  the 
coefficient  matrix  for  the  unknown  incremental  values  as  the  successive 
approximations  approach  the  true  eigenvalue  and  eigenvector. 

The  method  presented  in  this  chapter  overcomes  this  shortcoming.  To 
accomplish  this,  all  eigenvectors  corresponding  to  multiple  or  close  eigen¬ 
values  are  found  together.  As  in  the  method  of  Chapter  2,  this  method 
yields  the  eigenvalues  and  corresponding  eigenvectors  at  the  same  time. 

The  essence  of  the  method  consists  first  in  finding  the  subspace 
spanned  by  the  eigenvectors  corresponding  to  multiple  or  close  eigenvalues. 
The  subspace  is  found  using  the  Newton-Raphson  technique  in  a  way  suggested 
by  the  Robinson-Harris  method  [44],  If  the  eigenvalues  of  interest  are 
multiple,  any  set  of  independent  vectors  spanning  subspace  are  the  true 
eigenvectors,  but  if  the  eigenvalues  are  merely  close  together,  the 
vectors  must  be  rotated  in  the  subspace  to  find  the  true  eigenvectors.  The 
eigenvalues  are  obtained  as  a  by-product  of  the  process  of  finding  the  sub¬ 
space  and  any  subsequent  rotation.  In  this  method,  any  number  of  close 
eigenvalues  or  an  eigenvalue  of  any  multiplicity  can  be  found  together  with 
the  corresponding  eigenvectors. 


«<W4«1 
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The  theoretical  background  of  the  method  is  presented  in  Section  3.2 
The  iterative  scheme  for  finding  the  subspace  of  the  eigenvectors  corre¬ 
sponding  to  multiple  or  close  eigenvalues  is  given  in  Section  3.3.  The 
additional  treatment  required  for  close  eigenvalues  and  corresponding 
eigenvectors  is  the  subject  of  Section  3.4.  The  convergence  rate  and  the 
number  of  operations  per  iteration  are  given  in  Section  3.5. 

3.2  Theoretical  Background 

Let  us  consider  the  system  treated  in  Chapter  2,  i.e., 

Ax.j  =  Bx..  (i  =  1,  2,  .  .  ,  ,  n)  (3.1) 

where  A  and  B  are  symmetric  matrices  of  order  n,  and  B  is  positive  definite. 
The  x.  are  eigenvectors,  and  the  a^  eigenvalues  in  the  order  Aj<  a2<  . <a 

Let  a  set  S  consist  of  s  integers  p.  (j  =  1,  2,  .  .  .  ,  s),  that  is, 

w 

S  =  [p^,  P2»...,PS]  where  1  <  Pj<  n.  The  s-dimensional  subspace  spanned  by 
the  eigenvectors  x.  (jeS)  where  none  of  the  corresponding  eigenvalues 

J 

A.  (jeS)  are  close  or  equal  to  eigenvalues  A.  (i^S)  is  denoted  by  R.  Let 

J  * 

us  take  s  vectors  y .  (jeS)  which  are  orthonormal  with  respect  to  B  and  are 

J 

in  the  neighborhood  of  the  subspace  R.  This  means  that  if  the  vector  y.  is 

J 

expanded  in  a  series  of  true  eigenvectors  x^  (i  =  1,  2,...,n) 

n 

yj  =  L  xi  (jeS)  (3.2) 

i-1 


« 
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then,  the  following  relations  must  be  met: 


(jeS) 


(3.3) 


Hence,  a  vector  y.(jeS)  needs  not  be  close  to  one  of  the  x.(jeS). 

J  J 

With  the  above  definitions,  the  subspace  R  of  the  eigenvectors  x.(jeS) 

J 

is  characterized  by  the  following  constrained  stationary- value  problem:  find 
the  stationary  values  of 


subject  to 


n 

W  =  z  JjT  Ajj 

jeS 


(3.4) 


yl  By .  =  6 .  . 
i  J  1J 


(i,  jeS) 


(3.5) 


where  6..  is  the  Kronecker  delta,  i.e.,  6 . •  =  1  for  i  =  j,  and  6..  =  0  for 

I J  I  J  I J 

i  f  j.  The  function  w  could  be  regarded  as  a  sum  of  Rayleigh  quotients  of 
the  vectors  y.,  since  by  Eq.  (3.5)  the  denominators  of  the  Rayleigh  quotients 

J 

are  equal  to  unity.  The  important  result  that  the  stationary  property 

characterizes  the  subspace  R  is  proved  as  Theorem  1  of  Appendix  C. 

The  stationary- value  problem  may  be  treated  by  the  method  of  Lagrange 

multipliers.  Introducing  the  undetermined  multipliers  y . .  (i,jeS)  and  letting 

'  J 

y. .  =  y  ..  (see  Eq.  (3.5)),  we  have  the  Lagrangian 

I  J  J  I 


L  = 


*i  A*i  ' 


ieS 


ieS  jeS 


■J 


Wij  (yi  Byj  “  6ij) 


(3.6) 
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The  problem  of  Eqs.  (3.4)  and  (3.5)  is  equivalent  to  that  of  solving  the 
unconstrained  stationary-value  problem  for  the  Lagrangian  L.  The  problem 
is  solved  setting  the  first  partial  derivatives  of  L  with  respect  to  the 
unknowns  y.  and  u..  equal  to  zero,  i.e., 

v  *  J 


~  -  o  '•  Ay< 


■  E  ^  By,  (JeS) 


(3.7) 


-iL-  =0  ;  yl  By .  =  6 . . 


(i,  JeS) 


(3.8) 


Introducing  the  following  notation 


Y  =  [y  »  y  ,«».»  y  3 
P1  p2  ps 


aj =  {\r  yp2j,,"’Wp3j)  (j  =  ?v P2 . P$) 


D  =  (d  ,  d  »...»d  ) 

P1  p2  ps 


(3.9) 


we  can  write  Eq.  (3.7)  in  matrix  form  as 


Apj  =  BY3j  (j  =  Pj*  • *PS ) 


(3.10) 


or  collectively 


AY  =  BYO 


(3.11) 
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In  the  same  way,  Eq.  (3.8)  can  be  written  as 

YTBY  =  I  (3.12) 

where  I $  is  the  unit  matrix  of  order  s.  Hence,  the  subspace  R  of  the 
desired  eigenvectors  can  be  found  by  solving  Eqs.  (3.11)  and  (3.12).  Note 
that  Eqs.  (3.11)  and  (3.12)  are  nonlinear  in  D  and  Y  and  that  there  are 
s  (s  +  l)/2  scalar  unknown  elements  in  D,  since  D  is  symmetric,  and 
s  (s  +  l)/2  independent  equations  in  Eq.  (3.12).  In  the  next  section,  the 
solution  of  Eqs.  (3.11)  and  (3.12)  in  the  special  case  that  (jeS)  are  all 
multiple  or  close  eigenvalues  will  be  discussed. 

3.3  The  Iterative  Scheme 

In  this  section,  the  application  of  the  Newton- Raphs on  technique  to 
the  solution  of  Eqs.  (3.11)  and  (3.12)  for  multiple  or  close  eigenvalues 
and  their  corresponding  eigenvectors  will  be  presented.  To  simplify  the 
notation  in  this  discussion,  we  take  the  set  S  =  [l,2,...,s],  that  is,  the 
s  lowest  eigenvalues  are  close  together,  or  the  multiplicity  of  the  lowest 
eigenvalue  is  s.  It  should  be  emphasized  that  this  is  not  restrictive,  and 
the  procedure  is  perfectly  applicable  to  multiple  or  close  eigenvalues  in 
any  range. 

Assume  that  the  initial  values  for  D  and  Y,  and  Y^  are  available 
(the  solution  for  the  initial  values  will  be  discussed  in  Chapter  4). 
Furthermore,  we  assume  that  the  initial  vectors  in  Y^  are  in  the  neighbor¬ 
hood  of  the  subspace  of  the  eigenvectors  X  =  [x^ .x,,,. . . ,x$]  and  that  they 
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have  been  orthonormal ized  with  respect  to  the  matrix  B,  i.e.,  Y^BY^  =  I$. 
With  the  above  assumptions,  we  now  apply  the  Newton-Raphson  technique  to  the 
solution  of  Eqs.  (3.11)  and  (3.12).  For  the  general  kth  iteration  step,  let 


d  .(k+1)  =  d.^  +  Ad.^ 
J  J  J 


r.(k+D  =  zM 


(3.13) 


-  (k)  -  (k)  -  (k)  -  (k) 

where  Ad.'  '  and  Ay.  '  are  unknown  incremental  values  for  d.  '  and  y.  . 

J  J  J  J 

Introducing  Eq.  (3.13)  into  Eqs.  (3.10)  and  (3.12)  and  neglecting  the 

-  (k) 

nonlinear  terms,  we  obtain  the  linear  simultaneous  equations  for  Ad.'  '  and 

AAy,^  -  BY^Ad.^  =  -  Ay/k^  +  BY^d/k)  +  BAY^d.^  (3.14) 

J  J  J  J  J 


y(k)  By(k)  +  2Y(k)  gAY^k^  «  I 


(3.15) 


By  Theorem  3  of  Appendix  C,  if  the  (j  =  1,2 . s)  are  multiple  or 

J 

close  eigenvalues,  the  off-diagonal  elements  of  D  are  zero  or  very  small 

compared  with  its  diagonal  ones,  thus  the  last  term  of  Eq.  (3.14)  may  be 

-  (k) 

approximated  by  u..B  Ay.'  ' ,  yielding 

J  J  J 

(A  -  p..(k^B)  Ay.(k)  -  BY(k)Ad,(k)  =  -  Ay^k)  +  BY(k)d.(k)  (3.16) 

Jj  J  J  J  J 


y(k)  gy^k^  =  I 


(3.17) 


Let  us  take 
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Then,  Eq.  (3.15)  becomes 


Y^  BAy.^  =  0 
J 


(3.18) 


which  is  the  condition  that  the  incremental  vectors  be  orthogonal  to  the 
current  vectors  with  respect  to  B.  If  the  computational  scheme  is  slightly 
altered  so  that  the  latest  y.^  is  used  at  all  times,  the  orthogonality 
condition  is  satisfied  automatically  provided  that  the  initial  vectors  y. ^ 

are  orthogonal.  What  this  means  is  that  we  use  y.^  (i  =  1,2 . j  -  1) 

for  the  computation  of  y.^k+1^. 

J 

The  final  equations  to  solve  for  Ad.^  and  Ay.^  are  Eqs.  (3.16) 

3  3 

and  (3.18)  along  with  the  orthonormality  condition,  Eq.  (3.17).  These 
equations  can  be  written  in  matrix  form  as 


A  -  -  BY^ 

J  J 


y(k)  g 


4d<k> 

.  J  _ 

1 

o 

(3.19) 


where 


f  =  Ay . ^ k ^  -  BY^  d.^ 
3  3  3 


(3.20) 


The  coefficient  matrix  for  the  unknowns,  d.^  and  y.^,  is  symmetric. 

3  3 

Furthermore ,  It  is  nonsingular,  as  is  shown  in  Appendix  A,  Thus,  Eq.  (3.19) 
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-  (kl  -  m  -  fk+n 

can  be  solved  for  Ad.  '  and  Ay.v  ,  yielding  improved  values,  d.  '  and 

J  J  J 

yj(k+1)  from  Eq.  (3.13). 

The  algorithm  using  Eq.  (3.19)  requires  a  new  triangularization  in 
each  iteration,  since  the  coefficient  matrix  is  changed  in  each  iteration. 

It  therefore  seems  useful,  as  in  Chapter  2,  to  substitute  (A  -  y..^B)  for 

.  .  J  J 

( k ) 

(A  -  Pjj  'B)  in  Eq.  (3.19)  in  order  to  save  computational  effort  in  the 
solution.  That  is,  the  basic  equations  for  the  increments  are  taken  as 


A-^(0)‘ 

>- 

CO 

1 

■  ^ 

II 

-  (k)’ 
‘  rj 

-  B 

o 

.3  .<» 

m  J  m 

l - 

o 

( _ 

(3.21) 


where  the  residual  vector  r.^  is  defined  as  in  Eq.  (3.20).  The  coefficient 

J 

matrix  in  Eq.  (3.21)  is  also  symmetric  and  nonsingular  (Appendix  A).  The 
equation  (3.21)  was  obtained  discarding  a  small  linear  term  (u^.^  - 
v»i/0^)B  y of  Eq.  (3.19).  The  procedure  using  Eq.  (3.21)  requires  only 

J  J  J 

partial  triangularizations  in  each  iteration,  since  only  the  vectors  in 
(k) 

Y'  are  changed,  reducing  the  number  of  operations  per  iteration.  The  pro- 

cedure  depends,  for  its  convenience,  on  the  decoupling  of  the  Ay.'  '  for  the 
-  (k) 

s  vectors  y  .  (i=l,2,...  ,s).  The  decoupling  was  possible  only  because 

J 

the  small  linear  terms 
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(see  Eq.  (3.14))  could  be  dropped  for  X.  (j=l,  2,  .  .  .  ,  s)  all  close 
together.  Experience  with  Eq.  (3.21)  for  A.  (j=l ,  2,  .  .  .  ,  s)  which  are  not 

J 

close  together  indicates  that  satisfactory  results  cannot  be  obtained. 

Note  that  if  s  =  1,  Eqs.  (3.19)  and  (3.21)  are  equivalent  to  the 

equations  used  for  distinct  eigenvalues  and  corresponding  eigenvectors: 

Eq.  (3.19)  becomes  Eq.  (2.6),  the  equations  used  in  the  Robinson-Harris 

method,  and  Eq.  (3.21)  becomes  Eq.  (2.7),  used  in  the  proposed  method. 

(k)  -  (k) 

With  sufficient  large  k,  the  incremental  values  a3.v  '  and  Ay.v  ' 

J  J 

will  vanish.  Then,  from  Eq.  (3.21) 


lim  r.(k)  =  lim  (Ay,(k)  -  BY(k)  d-(k))  =  0 
k-*»  ^  k-*»  ^  ^ 


(3.22) 


Letting 


lima(k) 

dj  “  k~  ^ 


v  =1in,yi 
yj  k-*»  J 


(k) 


(3.23) 


we  write  Eqs.  (3.22)  and  (3.17)  as 


AY  =  BYD  (3.24) 

YTBY  =  I$  (3.25) 

where  Y  =  (yj,^*  •  •  •  »y$) »  an d  0  =  (dj,^,. . .  ,d$).  By  Theorem  3  of  Appendix  C, 

if  the  eigenvalues  Aj  (j=l,2 . s)  are  multiple,  the  values  of  the  off- 

diagonal  elements  of  D  are  all  zero,  and  its  diagonal  elements  have  an  equal 


B 
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value  which  is  the  desired  multiple  eigenvalue.  Moreover,  the  vectors  in 
Y  are  the  corresponding  eigenvectors.  However,  if  the  eigenvalues  are 
close  but  not  equal,  additional  operations  are  required  to  find  the  desired 
eigenvalues  and  eigenvectors.  These  additional  operations  are  the  subject 
of  the  next  section. 

3.4  Treatment  of  Close  Roots 

Once  the  converged  solution  D  and  Y  has  been  found  by  the  algorithm 
described  in  the  previous  section,  but  the  values  of  the  off-diagonal 
elements  of  D  are  not  zero,  the  vectors  in  Y  are  rotated  in  the  subspace 
of  Y  to  find  the  true  eigenvectors.  A  rotation  matrix  is  found  by  solving 
a  small  eigenvalue  problem.  Furthermore,  the  eigenvalues  of  the  small 
eigenvalue  problem  are  the  desired  eigenvalues.  The  derivation  of  the 
small  eigenvalue  problem  is  as  follows.  The  system  with  the  s  eigenvectors 
in  X  =  [xrx2,...,xs]  and  corresponding  eigenvalues  in  a  =  diag  (Xj.Xg,-..,: 
may  be  written  as 

AX  =  BXA  (3.26) 

where  A  and  B  are  symmetric  matrices  of  order  n.  Now,  let 

X  =  YZ  (3.27) 

where  Z  is  the  unknown  rotation  matrix  of  order  s.  Introducing  Eq.  (3.27) 
into  Eq.  (3.26),  we  get 


Z 
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Postmulti plying  Eq.  (3.24)  by  the  matrix  Z  yields 

AYZ  =  BYDZ  (3.29) 

Premultiplying  Eqs.  (3.28)  and  (3.29)  and  using  Y^BY  =  I$  of  Eq.  (3.25),  we 
obtain  the  special  eigenvalue  problem  of  order  s 

DZ  «  Za  (3.30) 

where  D  is  the  converged  solution  found  by  the  algorithm  of  the  previous 
section.  The  matrix  0  is  symmetric  (see  Eq.  (3.24))  and  of  order  s,  the 
number  of  close  eigenvalues,  which  is  usually  small.  The  absolute  values 
of  the  off-diagonal  elements  of  D  are  small  compared  with  those  of  its 
diagonal  elements  (see  Appendix  C).  The  eigenvalue  problem,  Eq.  (3.30)  can 
be  easily  solved  by  any  suitable  technique  such  as  Jacobi's  method  [31,51], 
yielding  the  desired  eigenvalues  in  A  (x^,  X,,,  .  .  .  ,  x$)  and  the  matrix  Z, 
which  in  turn  gives  the  eigenvectors  X  by  Eq.  (3.27).  The  number  of  oper¬ 
ations  required  for  the  solution  of  Eq.  (3.30)  is  very  small  compared  with 
that  of  Eq.  (3.21),  since  s  is  small. 

3.5  Convergence  Rate  and  Operation  Count 

In  this  section,  the  convergence  rates  of  a  multiple  eigenvalue  and 
the  corresponding  eigenvectors  found  in  Appendix  B  will  be  sumnarized.  For 
convenience,  we  assume  that  the  lowest  eigenvalues  are  multiple,  i.e. , 

if  _  /  l  \ 

x  =  Xj  s  x^  =  ...  *  x$.  Let  the  approximate  eigenvectors  y.  '  (j  =  1 ,  2, . . . ,s) 
be  expanded  in  terms  of  the  eigenvectors  x..  ( i  =  1 ,  2,  .  .  .  ,  n),  i.e.. 


i  =  l 


(k) 

where  .  is  a  scalar  representing  the  components  of  the  eigenvector  x.. 
on  y.^.  If  denotes  the  error  in  and  8.^  the  error  in  y.^, 

J  J  JJ  J  J 


then  they  may  be  defined  by 
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m 

i 

m 

t  can  be  seen  from  Eqs.  (3.34)  and  (3.35)  that  the  eigenvalues  and  the 
'*’*  corresponding  eigenvectors  converge  linearly.  However,  the  eigenvalues 

converge  much  more  rapidly  than  the  eigenvectors. 

»  * 

The  number  of  operations  Np  required  for  finding  multiple  or  close 
w ,  eigenvalues  and  the  corresponding  eigenvectors  is  calculated  in  Table  1. 

This  number  is 

*  * 

Np  =  pn  (m^  +  3mg  +  2m^  +  2)  +  Tpn  [(s  +  4)  mfl  +  2m|i  +  Jj-  (s^  +  7s  +  4)]  (3.37) 

where  s  is  the  multiplicity  of  an  eigenvalue  or  the  number  of  close  eigen- 
,  values,  and  Tp  is  the  total  number  of  iterations  required  for  a  solution. 

**  It  can  be  seen  that  if  s  =  1,  the  number  of  operations  is  equal  to  the  number 

of  operations  required  for  finding  a  simple  eigenvalue  and  the  corresponding 

*  * 

eigenvector  (see  Eq.  (2.14)). 

*• 
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4.  APPROXIMATE  STARTING  EIGENSOLUTION 


4.1  General 

The  iterative  methods  described  in  the  previous  chapters  begin  with  an 
approximate  starting  eigensolution.  In  this  chapter,  a  procedure  to  find 
the  starting  solution  is  presented.  The  approximate  starting  solution  of 
an  eigenvalue  problem  is  often  available  either  as  the  final  answer  in  some 
approximate  methods  or  as  an  intermediate  result  in  other  iterative  methods. 

Numerous  methods  for  approximate  solutions  have  been  developed.  These 
include  static  or  dynamic  condensation  [2,3,25,28,34,42],  Rayleigh-Ritz 
analysis  [48,51],  component  mode  analysis  [9,51],  and  related  methods  sum¬ 
marized  by  Uhrig  [50].  In  all  these  methods,  the  approximate  solution  is 
found  in  a  single  step,  and  not  in  an  iterative  process.  Hence,  automatic 
improvement  of  the  solution  is  not  built  into  the  procedure.  Moreover,  the 
success  of  the  methods  depends,  to  a  great  extent,  on  the  engineer's  judg¬ 
ment,  which  is  difficult  to  incorporate  into  an  automatic  computer  program. 

Another  possible  way  for  finding  the  approximate  solution  is  to  take 
the  intermediate  results  from  other  iterative  methods  such  as  a  method 
combining  the  Gram-Schmidt  orthogonal ization  process  [51]  with  simultaneous 
iteration  method  or  combining  Rayleigh-Ritz  analysis  [6,9,11,29,32,49]  with 
simultaneous  iteration  method.  The  latter  combined  method  is  sometimes 
called  the  "subspace  iteration  method"  [6,9].  The  subspace  iteration  method 
is  used  here  to  find  approximate  starting  solutions  because  it  has  a  better 
convergence  rate  than  most  others.  The  method  itself  turns  out  to  require 
selecting  starting  vectors.  However,  a  scheme  to  find  starting  vectors  for 

1 
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the  subspace  iteration  method  has  been  well  established  and  is  fairly  routine 
(see  Section  4.2.2).  In  the  next  section*  the  subspace  iteration  method  will 
be  discussed. 


4.2  Subspace  Iteration  Method 
4.2.1  The  Iterative  Scheme 

The  subspace  iteration  method  is  a  repeated  application  of  the 
classical  vector  iteration  method  (power  method)  and  Rayleigh-Ritz  analysis. 
Suppose  that  the  p  smallest  eigenvalues  (i  =  l,2,...,p)  and  corresponding 
eigenvectors  xi  are  required  and  that  we  have  p  initial  independent  vectors 
*.(0)  (i  =  1,2,... ,p)  spanning  a  p-dimensional  subspace  in  the  neighborhood 
of  the  subspace  of  the  desired  eigenvectors. 

If  the  approximate  eigenvectors  and  corresponding  eigenvalues  after  k 

iterations  are  denoted  by  x^k^  and  \.^k\  x2^ . x  ^], 

and  =  diag  (A^k^,  a2^  . .  ,Ap^) ,  the  subspace  iteration  method  for 

the  kth  iteration  may  be  described  as  follows: 

(i)  Find  the  improved  eigenvectors  =  [y^k\  y,/k^ ,. . .  ,y 
by.  the  simultaneous  inverse  iteration  method; 

AY(k)  =  bx^"1)  (4.1) 

(ii)  Compute  the  projections  of  the  operators  A  and  B  onto  the 

(kl 

subspace  spanned  by  the  p  vectors  in  Yv  ' ; 


J( k)  _  y(k)T  AY(k) 
g(k)  =  y(k)T  BY(k) 


(4.2) 


where  A  and  B'  are  pxp  symmetric  matrices. 

(iii)  Solve  the  eigenvalue  problem  of  reduced  order  p  for  the 

eigenvalues  in  =  diag  »•  ■  •  > Ap^  and  the 

eigenvectors  in  =  [z^k^,  . 5  ^]; 


J(k)  z(k)  =  g(k)  z(k)  D(k) 


(4.3) 


Then, 


(iv)  Find  an  improved  approximation  to  the  eigenvectors; 


x00  =  Y(k)  z(k) 


lim  Dv  '  =  diag  (Xj.  *2 . xp) 


(4.4) 


lim  X<k>  =  [x.,  x„...  ,x  ] 

k-x-  l  2  p 


(4.5) 


Note  that  Eqs.  (4.2)  through  (4.4)  represent  a  Rayleigh-Ritz  analysis  with 
the  vectors  in  as  the  Ritz  basis  vectors,  which  results  in  X^k\  the 
best  approximation  to  the  true  eigenvectors  in  the  subspace  of  Y^. 

More  rapid  convergence  can  be  obtained  by  taking  more  iteration  vectors 
than  the  number  of  eigensolutions  sought.  However,  the  more  starting 
vectors  are  taken,  the  more  computational  effort  is  required  per  iteration. 
As  an  optimal  number  of  iteration  vectors,  q,  q  =  min  (2p,  p  +  8)  has  been 
suggested  [6,9]. 
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To  find  eigenvalues  within  a  given  range  a  <  p  <  b  and  the  correspond¬ 
ing  eigenvectors,  we  may  use,  instead  of  Eq.  (4.1),  the  inverse  iteration 
with  a  shift  [32]: 

(A  -  uB)  Y*k)  =  BX(k_1)  (4.6) 

where  p  is  a  shift  and  can  be  taken  as  (a  +  b)/2.  It  is  clear  from  Eq.  (4.6) 
that  the  eigenvectors  corresponding  to  the  eigenvalues  in  the  vicinity  of  a 
shift  p  will  converge  rapidly.  However,  the  convergence  of  other  eigen¬ 
vectors  may  be  slower  than  when  the  shift  is  not  applied,  since  as  a  result 
of  the  application  of  the  shift,  the  absolute  values  of  some  shifted  eigen¬ 
values  may  become  closer. 

4.2.2  Starting  Vectors 

The  number  of  iterations  required  for  convergence  depends  on  how 
close  the  subspace  spanned  by  the  starting  vectors  is  to  the  exact  subspace. 
If  approximations  to  the  required  eigenvectors  are  already  available,  e.g. , 
from  a  previous  solution  to  a  similar  problem,  these  may  be  used  as  a  set 
of  starting  vectors.  If  not,  we  may  use  one  of  the  schemes  for  generating 
starting  vectors  which  have  been  proposed  as  effective  [6,11,32,47]. 

The  scheme  for  establishing  the  starting  vectors  proposed  by  Bathe  and 
Wilson  [6,9]  is  used  here  because  of  its  simplicity  and  effectiveness.  The 
scheme  may  be  described  as  follows.  The  first  column  of  BX^  in  Eq.  (4.1) 
is  formed  simply  from  the  diagonal  elements  of  B.  That  is,  if  BX^  is 
denoted  by  C, 

c.j1  ~  b^  (i  =  1,2,. . .  ,n)  (4.7) 


This  assures  that  all  mass  degrees-of- freedom  are  excited  in  order  not  to 
miss  a  mode  [6,9].  The  next  (q-1)  columns  in  C  may  each  have  all  zeros 
except  for  a  certain  coordinate  where  a  one  is  placed.  These  coordinates 
are  found  in  the  following  way.  First,  compute  the  ratios  a^/b.^ 

(i  =  1,2 . n)  and  take  the  (q-1)  s.'s  (j  =  1,2,..., q-1)  such  that  the 

J 

absolute  values  of  the  ratios  a^/b..  for  i  (i  =  Sj,  s2,...,sq_j)  are 
smallest  over  all  i.  Then, 

ci,  j-l  '  1  for  1  '  sj  (1  '  1‘2 . n) 

=0  for  i  f  s.  (j  =  1,2,. ...q-1)  (4.8) 

J 

If  the  absolute  values  of  the  ratios  are  close  or  equal,  then  it  was  recom¬ 
mended  [6,9]  that  the  s.'s  (j  =  1,2 . q-1)  be  chosen  so  that  they  are  well 

J 

spaced. 

4.2.3  Convergence  Rate,  Operation  Count,  and  Estimation  of  Errors 

With  an  adequate  choice  of  the  starting  vectors,  the  subspace 

iteration  method  gives  good  approximations  to  the  exact  eigenvalues  and 

eigenvectors  even  after  only  a  few  iterations.  However,  the  subsequent 

convergence  is  only  linear  with  the  rates  of  convergence  equal  to 

o 

Ai^Xq+l  ^  f°r  the  eigenvector  and  ( Ai/ xq+i )  f°r  the 

corresponding  eigenvalue.  These  ratios  indicate  that  for  the  higher  eigen¬ 
value  convergence  is  slower.  Hence,  the  convergence  of  the  pth  mode  control 
the  termination  of  the  iteration  process. 
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One  of  the  most  important  indicators  of  the  effectiveness  of  numerical 
methods  is  the  total  number  of  operations  required  for  finding  a  solution, 
which  depends  on  both  the  rate  of  convergence  and  the  number  of  operations 
per  Iteration.  This  number  for  the  subspace  iteration  method,  N$,  (see 
Table  1)  may  be  expressed  by 

N$  =  Tsqn  (2ma  +  4mb  +  2q  +  4)  +  n  (m^  +  3ma  +  mb  +  1) 

where  m„  and  m.  are  the  half  band-widths  of  A  and  B,  and  T  is  the  total 
number  of  iterations  required  for  the  solution. 

The  total  number  of  iterations  T  ,  depends  on  the  rate  of  convergence 
and  tolerances  of  the  errors  in  approximate  eigenvalues  and  eigenvectors. 
Bathe  and  Wilson  [6,9]  suggested  use  of  the  following  formula  for  the  esti¬ 
mation  of  errors  in  the  ith  eigenpair  at  the  kth  iteration: 


(4.10) 


-  (k) 

where  r^'  '  = 


=  (A  -  xi ^B)  x.^. 


The  error  estimated  by  Eq.  (4.10)  is  a  function  of  both  the  approximate 

eigenvalues  and  eigenvectors.  However,  it  may  be  more  reasonable  to  estimate 

the  errors  in  approximate  eigenvalues  and  eigenvectors  using  separate  formulas 
(k)  (kl 

as  follows:  let  '  and  e..'  '  be  the  errors  in  the  ith  approximate  eigen- 

(k) 

value  and  eigenvector.  Then  y^v  may  be  estimated  by 


,  (k+1)  .  (k) 

00  xi  ~  Xi 


(i  =  1,2,. . . , p ) 


(4.11) 


X  » 
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For  the  estimate  of  0^k^,  we  find  the  incremental  vectors  Ax^k^  from  the 
relations 


x  (k+1> 

xi 


00- 

11 


(k) 


AX. 


(k) 


Then, 


00TB4j.(k) 


=  0 


(4.12) 


6i 


(k) 


AX 


(k) 


BAXi 


(k) 


1/2 


(4.13) 


If  some  of  the  approximate  eigenvalues  x.  (i  =  p1§  p2 . p$)  are  equal  or 

very  close,  we  may  then  compute  Ax^k^  from  the  relations 


x  (k+1) 

xi 


£ 

i=P^ 


a-.(k)i.(k)  +  Ax  (k) 
ij  xj  AXi 


(k) 


BA^.(k) 


( j  =  p j iPg » • • • »Pg ) 


(k) 


Ps 

£ 

j=px 


2  (k)-  (k) 
a1j  Xj 


1/2 


(4.14) 


For  the  purpose  of  comparison  of  the  proposed  methods  of  Chapters  2  and  3 
with  the  subspace  iteration  method,  the  errors  were  computed  using  Eqs. 
(4.11)  to  (4.14). 
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4.3  Starting  Solution  for  the  Proposed  Method 


The  intermediate  results  from  the  subspace  iteration  are  used  as  the 
starting  solutions  for  the  proposed  method.  During  the  subspace  iterations, 
the  errors  in  approximate  eigenvalues  and  corresponding  eigenvectors  can  be 
estimated  by  the  scheme  described  in  Section  4.2.3.  Furthermore,  these 
errors  can  be  used  for  estimating  the  number  of  iterations  or  the  number  of 
operations  required  for  the  solution  by  both  the  subspace  iteration  method 
and  the  proposed  method.  Hence,  it  is  possible  to  estimate  the  optimal 
number  of  iterations  to  be  carried  out  by  the  subspace  iteration  method. 

This  optimal  number  of  iterations  is  usually  one  or  two. 

Let  A*  and  x*  (i  =  l,2,...,p)  be  the  intermediate  solutions  from  the 
subspace  iteration  method  after  the  optimal  number  of  iterations.  Then,  if 
the  x*  are  well  separated,  a*  and  x*  can  be  taken  as  the  starting  solutions 
for  the  method  of  Chapter  2,  A^0^  and  x.. ^ .  However,  if  some  of  them, 
e.g.,  a*  (i  =  pj.p^,. . . ,p$)  are  equal  or  very  close,  a*  and  x*  are  taken  as 
the  starting  solution  for  the  method  of  Chapter  3  as 

(0)  -* 

V  =  xi 


for  i  f  j  (i,j  =  pj,p2 . ps)  (4.15) 


-  W—, ■  JB— I  II  il»  -r**^*-*.*, 


It  should  be  noted  that  from  Eqs.  (4.3)  and  (4.4),  the  iteration 
vectors  in  the  subspace  iteration  method  are  always  orthogonal i zed  with 
respect  to  B.  Therefore,  orthogonal ization  is  not  required  for  the  first 
iteration  of  the  proposed  method. 


5.  NUMERICAL  RESULTS  AND  COMPARISONS 


5.1  General 

The  relative  efficiency  of  the  methods  developed  in  this  study  is 
illustrated  in  this  chapter  by  the  numerical  results  of  the  free  vibration 
analyses  of  the  following  example  problems: 

(a)  Ten-Story,  Ten-Bay  Plane  Frame 

(b)  Two-Hinged  Circular  Arch 

(c)  Simply  Supported  Plate. 

The  problems  were  formulated  using  a  stiffness  method  for  the  plane  frame 
problem,  a  finite  difference  method  for  the  arch  problem,  and  a  finite  element 
method  for  the  plate  problem.  No  attempt  has  been  made  to  present  the 
solutions  of  eigenvalue  problems  of  very  large  order,  although  the  proposed 
method  is  developed  for  them.  However,  some  trends  can  be  inferred  from  the 
example  problems  presented  here. 

The  first  two  problems,  with  distinct  eignevalues,  were  solved  by  the 
method  discussed  in  Chapter  2  and  the  third  one,  with  multiple  or  close 
eigenvalues,  by  the  method  of  Chapter  3.  The  above  problems  were  also  solved 
using  the  Robinson-Harris  method  [44]  and  the  subspace  iteration  method 
discussed  in  Chapter  4.  The  results  are  summarized  in  Tables  2  through  5.  The 
numerical  results  given  here  are  shown  to  be  consistent  with  the  convergence 
estimates  of  Appendix  B. 

For  each  method,  the  total  number  of  operations  required  for  finding  the 

desired  eigenvalues  and  eigenvectors  to  the  same  accuracy  was  found.  These 

-4 

are  presented  and  compared  in  Table  6.  Although  a  tolerance  of  10  on  the 
eigenvalues  and  eigenvectors  should  be  sufficient  for  normal  requirements,  it 
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was  taken  as  10  ®  for  the  purpose  of  comparisons  of  the  convergence  character¬ 
istics  of  the  methods.  - 

The  numerical  computations  of  the  above  problems  were  performed  on  the 
CDC  CYBER  175  system  of  the  Digital  Computer  Laboratory  of  the  University  of 
Illinois,  Urbana,  Illinois. 

5.2  Plane  Frame 

The  ten-story,  ten-bay  plane  shown  in  Fig.  2  was  taken  as  an  example 
problem  in  order  to  test  the  method  of  Chapter  2' for  problems  with  distinct 
eigenvalues.  The  problem  was  formulated  by  a  stiffness  method  in  which  the 
axial  deformations  of  the  members  are  considered,  but  the  shear  deformations 
neglected  [40].  The  frame  with  three  displacements  per  joint  has  a  total  of 
330  degrees  of  freedom.  The  mass  matrix  is  the  consistent  mass  matrix  [4,5] 
with  a  maximum  half-bandwidth  of  35,  equal  to  that  of  the  stiffness  matrix. 

The  four  smallest  eigenvalues  and  their  corresponding  eigenvectors 
were  computed  by  the  proposed  method,  by  the  Robinson -Harr is  method,  and  by 
the  subspace  iteration  method.  The  results  are  given  in  Table  2.  For  the 
subspace  iteration  method,  ten  starting  vectors  were  formed  by  the  technique 
suggested  by  Bathe  and  Wilson  (see  Section  4.2.2).  The  starting  approximate 
eigenvalues  and  eigenvectors  for  the  proposed  method  and  for  the  Robinson- 
Harris  method  were  established  by  performing  two  cycles  of  subspace  iteration. 
Table  2  shows  that  even  the  eigenvalues  calculated  by  two  subspace  iterations 
are  already  accurate  to  three  figures.  However,  the  eigenvectors  are  accurate 
to  only  one  or  two  figures.  In  addition,  the  convergence  of  eigenvectors  by 
the  subspace  iteration  method  is  so  slow,  as  discussed  in  Section  4.2.3,  that 
12  iterations  were  required  for  the  convergence  of  both  eigenvalues  and  eigen- 
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vectors  to  the  indicated  tolerance.  The  proposed  method  and  the  Robinson- 
Harris  method  required  only  two  iterations  for  the  convergence  of  eigenpairs 
except  for  that  of  the  fourth  mode,  which  required  four  iterations  by  the 
proposed  method  and  three  iterations  by  the  Robinson-Harris  method. 

The  total  number  of  operations  to  solve  for  all  the  desired  eigenpairs  by 

the  proposed  method  is  3.50x10^;  by  the  Robinson-Harris  method,  4.57x10^'  and 

6  1 

by  the  subspace  iteration  method,  9.27x10  .  Therefore,  the  Robinson-Harris 
method  required  1.31  times  as  many  operations  as  the  proposed  method  did,  and 
the  subspace  iteration  method  required  2.78  times  as  many  operations,  as  shown 
in  Table  6. 

5.3  Arch 

A  uniform  90  degree  circular  arch  simply  supported  at  both  ends  was 
analyzed  for  in-plane  vibration  behavior.  The  arch  has  the  radius  a  and  the 
thickness  h,  and  the  ratio  a/h  =  20.  Mel  in  and  Robinson  [36]  investigated  the 
free  vibration  behavior  of  such  an  arch  as  a  part  of  a  study  of  vibrations  of 
a  simply  supported  cylindrical  shell  using  a  finite  difference  method.  The 
arch  was  divided  into  12  uniform  segments  giving  22  degrees  of  freedom.  The 
maximum  half-bandwidth  of  the  stiffness  matrix  is  four  and  the  mass  matrix  is 
a  unit  diagonal  matrix. 

The  problem  was  analyzed  for  the  three  smallest  eigenvalues  and  their 
eigenvectors  by  the  proposed  method,  by  the  Robinson-Harris  method,  and  by  the 
subspace  iteration  method.  The  results  are  summarized  in  Table  3.  Five  radial 
displacements  were  taken  as  master  displacements  for  the  iteration  vectors  of 
the  subspace  iteration  method.  Starting  approximate  eigenpairs  for  the  proposed 
method  and  the  Robinson-Harris  method  were  established  by  carrying  out  just 
one  cycle  of  the  subspace  iteration. 
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The  comparison  of  the  total  number  of  operations  for  each  method  is 

3 

given  in  Table  6.  The  proposed  method  needed  8.87x10  operations,  the 

3 

Robinson-Harris  method  9.77X.0  operations,  and  the  subspace  iteration  method 

4 

1.76x10  operations.  Hence,  the  ratio  of  the  total  number  of  operations  by 
the  Robinson-Harris  method  to  that  by  the  proposed  method  is  1.10,  and  this 
ratio  for  the  subspace  iteration  method  is  1.98. 

5.4  Plate  Bending 

A  plate  simply  supported  on  all  edges  was  analyzed  in  order  to  test  the 
method  presented  in  Chapter  3,  for  the  solution  of  eigenvalue  problems  with 
multiple  or  close  eigenvalues.  The  plate  has  the  lengths  a  and  b,  and  the 
theckness  h.  Two  special  cases  were  considered;  an  aspect  ratio  b/a  of  1.00 
and  b/a  equal  to  1.01.  The  first  case  gives  multiple  roots,  while  the  second 
one  gives  close  roots.  The  problem  was  formulated  by  a  finite  element  method, 
in  which  the  plate  was  divided  into  16  elements.  Each  unrestarined  node  has 
a  deflection  and  two  rotational  displacements,  giving  a  total  of  39  degrees  of 
freedom.  The  mass  matrix  is  the  consistent  mass  matrix  [4,5]  with  a  maximum 
half-bandwidth  of  16,  equal  to  that  of  the  stiffness  matrix. 

The  four  smallest  eigenvalues  and  corresponding  eigenvectors  were  computed 
for  both  cases  by  the  proposed  method  and  by  the  subspace  iteration  method. 

The  results  are  summarized  in  Tables  4  and  5.  the  deflection  at  each  node  was 
taken  as  the  master  degrees  of  freedom,  giving  nine  iteration  vectors  for  the 
subspace  iteration  method.  Only  one  cycle  of  subspace  iteration  was  performed 
for  the  proposed  method.  The  multiple  eigenvalues  of  the  square  plate  and  the 
close  eigenvalues  of  the  rectangular  plate  were  isolated  by  the  method 
discussed  in  Chapter  3. 
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The  total  number  of  operations  by  the  proposed  method  for  both  cases  is 

5  5 

1.27x10  and  by  the  subspace  iteration  method*  2.20x10  ,  as  shown  in  Table  6. 

Hence,  the  subspace  iteration  method  needed  1.73  times  as  many  operations  as 

the  proposed  method  did. 

5.5  Comparison  between  the  Theoretical  Convergence  Rates  and  Numerical 
Results 

It  was  shown  in  the  previous  chapters  that  in  the  proposed  method,  the 

convergence  of  eigenvalues  is  much  faster  than  that  of  eigenvectors.  Hence, 

the  convergence  of  the  eigenvectors  governs  the  termination  of  process,  when 

the  tolerances  on  the  eigenvalues  and  eigenvectors  are  same.  Comparison 

between  the  theoretical  convergence  rates  and  numerical  results  was,  therefore, 

carried  out  only  for  the  eigenvectors.  Comparisons  between  the  proposed 

method  and  subspace  iteration  method  are  given  in  Tables  7,  8,  and  9. 

(k+1)  (k) 

The  numerical  convergence  rates  were  computed  by  el  '/el  ,  where 
(k) 

el  '  is  the  error  on  the  ith  approximate  eigenvector  at  the  kth  iteration. 

These  errors  are  given  in  Tables  2  through  5,  showing  that  the  numerical 
convergence  rates  for  the  proposed  method  and  the  subspace  iteration  method 
increase  monotonically  to  approach  the  theoretical  convergence  rates  as  the 
number  of  iterations  increases.  A  typical  example  for  this  is  the  convergence 
rates  of  the  fourth  eigenvector  of  the  frame  problem,  as  shown  in  Table  7. 

The  number  of  iterations  for  this  mode  is  large  enough  to  provide  a  good 
comparison  between  the  theoretical  and  numerical  convergence  rates. 

Tables  7,  8,  and  9  show  that  in  the  proposed  method,  eigenpairs  converge 
much  faster  than  in  the  subspace  iteration  method.  Note  also  that  in  Table  9, 
the  numerical  convergence  rates  for  the  proposed  method  are  almost  same  as 
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those  rates  for  the  problem  with  double  roots.  Hence,  the  expressions  for  the 
theoretical  convergence  rates  for  multiple  eigenvalues  also  seem  applicable  to 
the  case  of  close  eigenvalues. 
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6.  SUMMARY  AND  CONCLUSIONS 

6.1  Summary  of  the  Proposed  Method 

Two  iterative  procedures  for  the  solution  of  linear  eigenvalue  problems 
for  systems  with  a  finite  number  of  degrees  of  freedom  were  discussed  in 
Chapters  2  and  3.  Chapter  2  developed  a  procedure  for  finding  distinct 
eigenvalues  and  the  corresponding  eigenvectors,  and  Chapter  3  dealt  with 
multiple  or  close  eigenvalues  and  the  corresponding  eigenvectors. 

For  distinct  eigenvalues  and  the  corresponding  eigenvectors,  the  Robinson- 
Harris  method  [44]  was  modified  to  save  overall  computational  effort  by  the 
use  of  a  "modified"  form  of  the  Newton-Raphson  technique.  The  modified  method 
reduces  both  the  number  of  operations  per  iteration  and  the  convergence  rates. 
However,  the  reduction  of  the  number  of  operations  generally  compensates  for 
the  disadvantage  of  the  decrease  of  the  convergence  rate,  reducing  the  total 
number  of  operations. 

The  procedure  in  Chapter  2  for  finding  a  distinct  eigenvalue  and  the 
corresponding  eigenvector  fails  if  the  eigenvalue  is  one  of  multiple  or  close 
eigenvalues,  because  the  matrix  involved  in  the  computation  become  ill-condi¬ 
tioned.  This  difficulty  has  been  overcome  by  the  new  method  of  Chapter  3.  In 
this  mehtod,  all  eigenvalues  close  to  an  eigenvalue  or  a  multiple  eigenvalue 
and  the  corresponding  eigenvectors  are  found  in  a  group.  In  other  words,  a 
subspace  spanned  by  the  approximate  eigenvectors  is  projected  by  iterations 
onto  the  subspace  of  the  exact  eigenvectors.  If  the  eigenvalues  are  multiple, 
the  vectors  spanning  the  subspace  are  exact  eigenvectors.  However,  if  the 
eigenvalues  are  close,  the  exact  eigenvectors  are  found  by  a  simple  rotation 
of  the  vectors  in  the  subspace.  The  rotation  matrix  is  found  from  a  special 
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eigenvalue  problem  of  small  order  s,  the  number  of  the  close  eigenvalues. 

The  eigenvalues  of  the  small  eigenvalue  problem  are  exact  eigenvalues  of  the 
original  system. 

The  above  procedures  of  the  successive  approximations  require  initial 
approximations  to  the  eigenvalues  and  eigenvectors.  These  are  available 
either  as  the  final  solution  in  some  approximate  methods  such  as  static  or 
dynamic  condensation  or  as  an  intermediate  result  in  an  iterative  method  as 
the  subspace  iteration  method  described  in  Chapter  4. 

6.2  Conclusions 

The  method  presented  in  this  study  is  very  efficient  for  finding  a  limited 
number  of  soltutions  of  eigenvalue  problems  of  large  order  arising  from  the 
linear  dynamic  analysis  of  structures.  The  features  of  the  method  are  summarized 
as  follows. 

(a)  The  method  ha;,  very  high  convergence  rates  for  eigenvalues 
and  eigenvectors.  The  method  is  more  economical  than  the 
subspace  iteration  method,  the  advantage  being  greater 

in  larger  problems.  For  comparable  accuracy,  a  ten-story 
ten-bay  frame  required  only  36 %  of  the  number  of  operations 
need  in  applying  subspace  iterations. 

(b)  A  transformation  to  the  special  eigenvalue  problem  is  not 
required.  Thus,  the  characteristics  of  the  given  matrices 
such  as  the  sparseness,  bandness,  and  symmetry  are  preserved, 
minimizing  the  storage  requirements  and  the  number  of 
operations. 
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(c)  Any  number  of  multiple  or  close  eigenvalues  and  their 
eigenvectors  can  be  found.  The  existence  of  the  multiple 
or  close  eigenvalues  can  be  detected  during  the  iterations 
by  the  method  of  Chapter  2. 

(d)  The  eigenvalues  in  any  range  of  interest  and  their 
eigenvectors  can  be  found,  if  approximations  to  the 
solution  are  known. 

(e)  The  solution  can  be  checked  to  determine  if  some  eigenvalues 
and  corresponding  eigenvectors  of  interest  have  been 
missed,  without  extra  operations. 

6.3  Recommendations  for  Further  Study 

Several  possible  areas  of  further  study  to  improve  the  proposed  method 
may  be  suggested. 

(a)  The  convergence  rate  may  be  improved  by  other  modifica¬ 
tions  of  the  successive  approximation  method  used  for 
the  proposed  method. 

(b)  Further  improvements  may  be  possible  for  the  method  of 
finding  an  initial  approximation  to  the  eigensolution, 
and  for  isolating  the  eigenvalues  and  their  eigenvectors 
which  may  be  missed  by  the  proposed  method. 

(c)  The  proposed  method  may  be  applied  to  other  practical 
problems  of  our  interest  such  as  a  stability  analysis  of 
structures. 

(d)  The  proposed  method  could  be  easily  extended  to  the  contin¬ 
uous  eigenvalue  problems  if  there  were  better  ways  of 
direct  estimation  of  their  eigensolutions. 
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APPENDIX  A 


NONSINGULARITY  OF  THE  COEFFICIENT  MATRICES 
OF  THE  BASIC  EQUATIONS 


Consider  the  basic  equations  (3.19)  used  for  finding  multiple  or  close 
eigenvalues  and  corresponding  eigenvectors  of  the  system 

Ax  =  x  Bx  (A. 1 ) 

in  which  A  and  B  are  symmetric  and  both  of  order  n,  and  B  is  positive 
definite. 

Let  the  coefficient  matrix  of  Eq.  (3.19)  be  denoted  by  F,  that  is 


F  = 


A  -  Ulj(0)B 


-  y^)  b 


BY 


(k) 


i  =  m,m+  1 . m+  s-1  <n 


(A. 2) 


where  (i  =  m,m+  1,. . .  ,m+  s-1)  are  initial  approximate  values  of  the 

multiple  or  close  eigenvalues  A.,  (i  =  m,m+ 1,. . .  ,m+ s-1) ,  and  the  s  vectors 
in  Y^  =  [ym^,  yITTf/^»-'*»ynifS_i^3  are  approximate  values  of  the  eigen¬ 
vectors  in  X  =  [xm  Note  that  F  is  an  (n+s)x(n+s) 

symmetric  matrix. 

The  determinant  of  F  is  a  continuous  function  of  the  approximate  eigen¬ 
value  and  eigenvectors.  Hence,  if  F  is  nonsingular  when  the  approximate 


1 _ *  * 
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values  in  F  become  the  exact  ones,  then  it  will  be  also  nonsingular  for 
close  enough  approximations.  It  is  therefore  sufficient  for  our  purpose  to 
use  the  exact  eigenvalue  and  eigenvectors  in  Eq.  (A. 2)  to  prove  non¬ 
singularity.  Let  us  take  m  =  1  for  the  convenience  of  the  following  presen¬ 
tation,  then  the  resulting  matrix  G  will  be 


A  -  x.B 


-xtb 


(i  -  1,2,.. ,,s) 


(A.  3) 


where  X  =  [Xj,x2 . x$]. 

To  find  the  determinant  of  G,  we  follow  the  idea  that  Robinson  and 
Harris  [44]  used  for  showing  the  nonsingularity  of  the  coefficient  matrix 
of  Eq.  (2.10),  that  is,  we  consider  the  eigenvalues  y's  and  corresponding 
eigenvectors  u's  of  the  system 


*_ 

Gu  =  yBu 


(A. 4) 


or  collectively 


GU  =  BUD 


(A.  5) 


where 


U  =  ( ux  *u2 *  *  "  ,Us 


iff 


C  I 
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D  =  d  i  d  cj  (X^,X2*.**»X^) 

I  =  unit  matrix  of  order  s 
s 

It  may  be  verified  by  direct  substitution  that  the  (n+s)  eigenvectors 
u's  and  their  corresponding  eigenvalues  y's  of  Eq.  (A. 4)  are: 

j  =  1,2,...  ,s 
k  =  s+l,s+2,...,n 


Y  :  <-  iTT  >  <ejj>  -  V  (A-6> 

J  J 

where  xi  and  x.  are  the  eigenvalues  and  eigenvectors  of  the  system  Ax  =  xBx. 

_*  * 

The  vectors  ej  and  e^  form  the  diagonal  matrices  E  and  E  such  that 


E  [o i >^2 » • • • *e^ ]  -  diag  (e>| i *^22* *  *  * *^ss^ 


*  _  ★  .*  -  i  «*|  .]  . 

E  -  [e1,e9,...,e  J  =  diag  (5 — ,  ....  « —  ) 
12s  en  e22  ess 


e . . 
JJ 


a.n  dA 


+  4 


j  1,2,. • .  ,s 


A..  =  A.  -  A. 
JT  J  1 


i>J  —  1  )2  y  •  <  •  )S 


(A.  7) 


-4 


i 
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From  Eqs.  (A. 5),  (A. 6),  and  (A. 7) 


det  G  =  (det  B*)  (det  D) 

=  (-l)s  (det  B)  n  (a.  -  a.)  (A. 8) 

k=s+l  *  1 

In  a  similar  way,  the  determinant  of  G  for  general  m  >  1  is 

*  n 

det  G  =  (-l)s  (det  B)  n  (a.  -  A.)  (A.9) 

k=l  K  1 
MS 

where  the  set  S  =  [m,m+l . m+s-1].  The  matrix  B  is  positive  definite, 

which  implies  that  det  B  >  0.  Thus,  if  A.  (ieS)  is  not  equal  or  close  to 
any  of  A^  (k  =  l,2,...,n;  MS),  the  determinant  of  G  is  never  equal  to  zero 
or  close  to  zero,  independently  of  whether  xi  (ieS)  are  close,  multiple,  or 
distinct. 

From  Eq.  (A. 2),  if  s  =  1,  the  matrix  F  becomes  the  coefficient  matrix 
of  Eq.  (2.10),  and  by  Eq.  (A.9),  the  determinant  of  the  matrix  can  be 
approximated  by 


n 

F  =  (-1)  (det  B)  n  (a.  -  A  )  (A. 10) 

k=l  K  m 
Mm 


Therefore,  if  a  ^  A__.  and  Am  f  A,..,  the  matrix  F  is  also  nonsingular, 
m  m-i  m  m+i 
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APPENDIX  B 

CONVERGENCE  ANALYSIS 


The  convergence  analysis  of  the  methods  introduced  in  Chapters  2  and  3 
will  be  presented.  The  eigenvalue  problem  we  deal  with  here  is 


Ax..  =  A.  Bx.j  (i  =  1,2,. ...n)  (B.l) 

in  which  A  and  B  are  symmetric  matrices  of  order  n,  and  B  is  positive 
definite.  The  eigenvectors  x..  (i  ■  l,2,...,n)  are  assumed  to  be  ortho- 
normalized  with  respect  to  B. 

B.l  Case  of  a  Distinct  Root 

Let  us  rewrite  Eq.  (2.7)  used  for  improving  approximate  values  of  a 
distinct  eigenvalue  A.  and  the  corresponding  eigenvector  x.  of  the  system 

J  J 

represented  by  Eq.  (B.l): 


(A  -  a.^B)  Ax.(k)  -  AA Bx.^  =  -  (A  -  A.^B)  x.(k) 
J  J  J  J  J  J 


(B.2) 


X.(k)  BAX.*k)  =  0 

J  J 


(B.3) 


(k)  -  (k) 

where  A.v  and  x.'  are  approximate  values  of  A.  and  x.  after  k  iterations 

v  v  J  J 

and  aa . ^ k^  and  ax.^  are  unknown  incremental  values  of  A.^  and  x/k^. 

J  J  J  J 

_  (|<\  _  (L) 

Let  the  approximate  eigenvector  x.'  and  the  incremental  vector  ax- 

J  J 


be  expanded  in  a  series  of  the  true  eigenvectors,  i.e., 
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x 


j 


00  . 


n 


i=l 


n 


)  AC  •  • 
^  1J 


i=l 


(B.4) 


in  which  c..^  and  ac..^  are  scalar  coefficients.  Since  the  x.^  is  in 
i  j  i  j  J 

the  vicinity  of  x.» 

J 


max 

i 

i/j 


00 

T*7 


e  <<  1 


(B.  5) 


The  errors  in  A. 

J 


00  and  i.(k),  y . ^ k^  and  e.^O,  may  be  defined  by 
J  J  3 


oo  « 


A.  -  A. 

_J _ 1 


00 


Xj 


«  1 


V 


1/2 


«  1 


(B.6) 


where  the  values  of  y.00  and  e.00  are  very  small  compared  with  unity.  If 

J  3 

the  vectors  c^*0  and  c.  k  are  defined  by 

J  J 
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-(k)'  _  (c  (k)  c  (k)  (kh 

cj  “  lclj  ’  c2j  *  * "  *cnj  1 


c*(k)  =  (0,  0,  Cjj(k),  0 . 0) 


(B.  7) 


_  (k) 

The  e.v  '  then  represents  very  closely  the  angle  between  the  vectors  c. 

J  J 


-*(k) 


(k+D  «  (k+l )» 


and  c.v  '  (see  Fig.  1).  The  task  here  is  to  estimate  and  e 

J  J  1 

the  errors  in  A.^k_1^  and  x.^k_1^. 

J  J 


-  T 


Let  us  substitute  Eq.  (B.4)  into  Eqs.  (B.2)  and  premultiply  by  x.  to 


obtain 


(>  -  i  <0>1  AC  <k>  -  A»  M  -  -  U  -  J  <k))  C  <k) 

Ui  xj  >  4cij  cij  Ui  ’  Sj 


(i  =  1»  2,. . . ,n) 


(B.8) 


Substitution  of  Eq.  (B.4)  into  Eq.  (B.3)  and  use  of  the  orthonormality 
of  the  eigenvectors  with  respect  to  B  results  in 


I  CyW  -  0 

i-1 


(B.9) 


(k)  (kl 

The  unknown  quantities,  AA .  '  and  ac.  .'  '  will  be  found  from  Eqs.  (B.8) 

sJ  '  J 


and  (B.9). 
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where 


h  =  max 
i 

7j 


illii 

Ai  -  Aj 


(o) 


(0T 


«  1 


3  =  max 
i 

i7j 


VLii 

Ai  "  Aj 


00 

m 


-  1  (i  =  1,  2,. . .  ,n) 


/L] 

Therefore,  from  Eq.  (B.ll),  the  AX.V  '  may  be  approximated  by 

J 


AA.^  =  x.  -  x<k>  +  (x.  -  x.^)e 

J  J  J  J  J 


or 


X <k+1>  *  X  <k>  ♦  AX  <k> 
J  J  J 


■  xj  +  (xj  -  xj(0)>6 


Substitute  Eq.  (B.15)  into  Eq.  (B.IO)  to  obtain 


-»(M  ■  -  sj<w  ♦  iiJLV-  (I  +  6)  c”(k) 

Ai  "  Aj 


'ij 


,(°) 

W 


'ij 


(B.14) 


(B.15) 


(B. 16) 


(i  -  1,  2,. . . ,n) 


(B.17) 
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c  .(k+1>  =  c  {k)  +  AC  -(k) 
Cij  cij  +  Acij 


X.  -  X  v/ 

=  J (1  +  e)  cfj(k)  (i  =  1,  2 . n)  (B.  18) 

xi  '  Xj 


from  which,  it  follows  that 


cjj(kn)  =  <» *  «>  cij<k) 


(B.  19) 


The  measure  of  the  error  in  x^k+1^,  e.^k+1^,  can  now  be  found: 

J  J 


n  ,  (k+1)  ,2  1/2 

(k+l)  .  r  hi _  \ 


^  7WJ 

i=l  '  jj 
ij'j 


X..-A<°V  /  cj»\2 


xi  “  xj 


1.1 

T^W 


<  h  0.' 


(B.20) 


where  h  is  given  in  Eq.  (B.14)  and  is  very  small  compared  with  unity.  To 
find  Y^k+1^,  the  measure  of  the  error  of  x^k+1\  we  use  Eqs.  (B.6)  and 
(B.16),  giving 


(k+1)  = 


X,  -  A. 


(k+1) 


A.  -  A. 


A .  Aj 


g 

J  J 


(B.21) 


by  which 


(k+2)  < 


A.  -  A. 


s2(k+l) 


(B. 22) 


Substitution  of  Eq.  (B.20)  into  Eq.  (B.22)  and  use  of  Eq.  (B.21)  results  in 


V.<k+2)  -  h2  y . (k+1) 

J  J 


(B. 23) 


Hence,  it  can  be  seen  from  Eqs.  (B.20)  and  (B.23)  that  the  jtn  eigenvector 

and  eigenvalue  converge  linearly  with  errors  multiplied  by  h  (h  <<  1  )  and 
2 

h  respectively  in  each  iteration. 

B.2  Case  of  a  Multiple  Root 

The  convergence  analysis  of  the  method  for  finding  a  multiple  eigen¬ 
value  and  the  corresponding  eigenvectors  of  the  system  given  in  Eq.  (B.l) 
will  now  be  presented. 


'»*******«*•' r« 
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For  convenience,  but  without  loss  of  generality,  the  s  lowest  eigen¬ 
values  are  assumed  to  be  equal,  and  the  eigenvalue  of  multiplicity  s  is 
denoted  by  ,  1 .e. ,  \*  =  Xj  =  \2  =  •  •  •  =  xs*  Let  us  rewrite  Eqs.  (3.21) 
and  (3.17),  which  are  the  basic  equations  for  improving  approximate  values 
of  the  multiple  eigenvalue  and  corresponding  eigenvectors,  i.e.. 


and 


(A  -  y..(0)B)  Ay.^  -  BY^  Ad.^  =  BY^  d/k^  -  A  y.^ 

JJ  J  J  J  J 


Y(k)  B  Ay.OO  =  0 
J 


(j  =  1,  2 . s) 


(j  =  1,  2,. . . ,s) 


(B. 24) 

(B. 25) 


y(k)T  BY(k)  =  J 


(B.26) 


where  I$  is  the  unit  matrix  of  order  s,  and 


Y(k)  =  til00.  ;2(k) . yslk)] 


d  <k>  =  (u  <k>  u  (k)  „  (k>) 

j  ulj  *  u2j  ’•‘•’wsj  > 


.3  (k)T  _  /A  (k)  .  (k)  A  (k) , 

Adj  -  (avjj  ,  Au2j  . AwSJ.  ) 


(B.27) 


( k) 

The  y..'  '  (j  =  1,  2 . s)  are  approximations  to  the  multiple  eigenvalue 

.  .  /  k  ) 

a  =  xi  3  x2  3  '  ’  •  =  xs»  and  the  '  (j  =  1»  2 . s)  are  approximations 
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to  the  true  eigenvectors  x.  (j  s  1,  2,...,s).  The  Ay.^  and  Ad.^  are 

-  (k)  -  (k)  J 

unknown  incremental  vectors  for  y.'  and  d.  . 

J  J 

-  fkl 

Let  the  approximate  eigenvectors  y.'  '  and  the  incremental  vectors 

_  (L\ 

Ay.  be  expanded  in  a  series  of  the  true  eigenvectors  x.  (i  =  1,  2 . n), 

J  ■ 

as  in  Eq.  (B.4),  i.e. , 


(k)  _ 


L  Cid 
i=l 


n 

Ay  j ( k }  =  Z  Acij(k)  (J  -  1.  2 . s)  (B.28) 

i=l 

(kl  -  (kl  Ik)  fkl 

Denoting  the  errors  in  u . .  and  y.'  '  by  y.  and  e.v  ,  we  have 

J  J  J  J  J 


where 


1 


9j 


e  (k) 
(k)  -  Krr 
“j 


«  1 


,  %  (  s  1  1/2 

.00  J  y  c2  (k)  V 


L 

li-l 


ij 


(k) 

j 


( i=s+l 


1 1/2 

,2  (k)  \ 

,j  J 


(B.29) 


(B. 30) 
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where 


A  ~  di  eg  (Xj»  ^2 i • • • »^p) 


I  =  unit  matrix  of  order  n. 
n 


Let  us  find  Ac,^  and  Ad;^  using  Eqs.  (B.35)  -  (B.37).  From  Eq.  (B.35) 

J  J 

AC.(k)  =  (A  -  „  <°>  I  J"1  (c(k)  d.(k+1)  -  A  l.W)  (B.  38) 

J  J  J  ''  si 

Substitution  of  Eq.  (B.38)  into  Eq.  (B. 36)  leads  to 


F(k)  a(k+l)=  -  (k) 


(B. 39) 


where 


F(k)  ,  c(k)T  ^  _  (0)  jj-1  c(k) 


JJ  n' 


g  (k)  =  c(k)T  (A  _  (0)  J  )-lA  -  (k) 


JJ  n* 


(B.40) 


( k) 

Note  that  Fv  '  is  a  symmetric  matrix  of  order  s.  Using  Eq.  (B.33),  we  can 
show  that 


F(k)  =  - - 1  •  ^  R^k)  (Is  +  E(k))  R^k) 

X  "  pjj 


(B.41) 


where 


R(k)  =  diag  (aj(k\  a2(k) . c.$(k)) 


(B.42) 
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and  the  (£,  m)t(l  element  of  the  symmetric  matrix  E^,  is 


(k)  _ 


mm 


-  aff 

;  _  .  to) 


<C.  <k> 
lm 


i=s+l  i  Hjj 


\am 


W 


(m  "  1»  2,. . .  ,s) 


xm 


(k)  . 


X*-  X, 


,'c  <k>  \  '  c.  \ 

1  i  im 


i=s+l  Ai  "  yjj 


tot  orj 


\  at' 


u 


m 


(k) 

wrj 


U,  m  =  1,  2,. . .  ,s;  it  m)  (B.43) 


The  (i  =  1,  2,. . .  ,s)  of  Eq.  (B.42)  are  defined  in  Eq.  (B.30).  By 


Eq.  (B.29) ,  the  absolute  values  of  (l,  m  -  1,  2,...,s)  are  very  small 


compared  with  unity,  thus 


r-(k)  ,  1  ...  ,  2(k)  2(k)  2(kh  ,D  nA^ 

F'  =  (— - diag  (oy4  ,  a24  ',...,os  )  ( B . 44 ) 


1  •  ujj 


00 


Similarly,  the  values  of  the  s  components  of  the  vector  g.  can  be  found 
from  Eq.  (B.40) ,  i .e. , 


2(k) 


JJ 


9ij<k)  ■  V  ~  ajTo)  -jj'0’  '>ij(k)  -  »•  2 . .  (  * J)  <8'45) 

A  "  Wjj 


rrssss 
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where 


(k)  _ 


JJ 


”  ^ii  U  \«4 


m=s+l  m  Hjj 


i  c 

w  Hferl 


i  j 


(k)..r  »  -  »*  IcJVMcJU' 


L  _ EL 


,  mi  i  f  mi 

W  7TkT~ I  (k r 


«+1  K'"'  MV 


(B.46) 


Since  d  (k+1)T  -  t  <k+1)  (k+1)  (k+1)  k  *•  .*•  / 

iince  ck  -  ,  p^j  PSJ-  by  definition  (see 

Eq.  (B.2  )),  Eqs.  (B.39),  (B.4  ),  and  (B.45)  result  in 


„  (k+1)  „  .*  .  (Oh  (k)  , .  ,  ,  . 

“jj  -»+(»-  UJJ  )  njJv  (0  =  1.  2 . s) 


l,ij<IC+1)  "  “oo<0>  "io<k)  <1>j  =1-  2 . Si  4  *  3) 

Substitution  of  Eq.  (B.43)  into  Eq.  (B.38)  results  in 


(B - 47 ) 


ACij 


(k)  =  1  '7  r  (k)  (k+1) 

7 ” ToT  Cim  pmj 
Ai  wjj  m=l 


a. 

i 


Ai  •  Ujj 


Wcij 


(k) 


X*  -  A 


Ai  "  PJj 


i _  c  (k)  + 

W)  cij  + 


A*  S  (°) 

_  r~’  ffi  n  (k)  I  Pmj  (kh 

A  _  „  (0)  u  t6mj  njj  *  nmj  ] 

Ai  ujj  m=l  A 


c  (k)  +  L-'-^  [1  +  0  (e2)]  c  W 

Xi  '  “jj 


(i  -  1 1  2,. . .  ,n) 


(B.48) 
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or 


c. . 
U 


,*  (0) 

(k«)  x  '  1,jJ  ,  c  M 

ToTcij 


(i  =  1 ,  2,. . .  ,n) 


Ai  •  "jj 


from  which,  since  A*  =  -  ■■■  ~  An>  it  follows  that 

c,  <ktl)  .  e.  <»  (1-1.2 . s) 

'  J  *  J 

Thus, 


9j 


(k+1)  _ 


£  (c.(k+1> 

i=s+l  J 


■ -,r!)2 

i  =  l 


1/2 


(B. 49) 


(B. 50) 


n  * 

A  -  p . 
r-  J 

(0) 

J, 

ivn’l 

,i=s+l  xi  "  uj 

(oj 

j 

1^/  J 

1/2 


(k)  « 


2  h  e.(k) 
J 


where  y.'  is  defined  in  Eq.  (B.30),  and 


h  =  max 
i 


x*  -  ^(0) 


(B. 51) 


(i  =  s+1,  s+2,.. .  ,n)  (B.52) 


To  find  yJk+^f  the  measure  of  the  error  of  p..^k+^»  Eqs.  (B.29) 

J  J 

and  (B.47)  are  used,  which  results  in 


X*  -  U  <k+1) 

( k+1 )  _  A  UJJ 


X*  -  u  (0) 

_ _ -  n  (k) 

x*  JJ 


(B. 53) 


where  is  given  in  Eq.  (B.46).  Its  absolute  value  is 

J  J 


(k)  <  Ai  _2(k) 

njj'  -  mi?x  - Toy  ej 

JJ  l  x.  -  yi1'  '  J 

s  J  J 


(i  =  s+1,  s+2,...,n)  (B.54) 


Therefore,  from  Eq.  (B.53) 


where 


VM)  -  i  e?(k) 


x*  -  ,  <« 


(B.55) 


5  =  max 


(i  =  s+1,  s+2,...,n)  (B.56) 


1  Ai  -  "jj 


From  Eqs.  (B.55)  and  (B.51), 


Y  (k+2)  _  fl2(k+l) 
Yj  "  C  9j 


h2  e2(k) 


h2  y.(k+1)  (j  =  1,  2 . s)  (B. 57) 

J 
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Hence,  it  can  be  seen  from  Eqs.  (B.5n)  and  (B.51)  that  the  multiple  eigen¬ 
value  and  the  corresponding  eigenvectors  converge  linearly  with  errors  multi- 
2 

plied  by  h  (h  <<  1)  and  h  respectively  in  each  iteration. 

From  Eqs.  (B.47)  and  (B.49), 


lim  =  A*  for  i  =  j  (i,j  =  1,  2,...,s)  (B.58) 

k^>  1J 

=  0  for  i  i  j 
and 

lim  c..<k>  =0  (i  =  S+l,  s+2,...,n; 

k-*»  1J 

j  =  1,  2,...,s)  (8.59) 


which  shows  that  as  k->®,  the  vectors  x.^  (j  =  1,  2,...,s)  span  the  subspace 

J 

of  x.  (j  =  1,  2,...,  s)  whose  corresponding  eigenvalue  is  multiple.  Thus, 

J 

the  x.(k)  (j  -  1,  2,...,s)  themselves  are  a  set  of  true  eigenvectors,  ortho- 

J 

normalized  with  respect  to  B. 


APPENDIX  C 


THE  BASIC  THEOREMS  ON  THE  CONSTRAINED  STATIONARY- VALUE  PROBLEM 

Three  theorems  used  for  the  development  of  the  method  for  finding 
multiple  or  close  eigenvalues  and  the  corresponding  eigenvectors  will  be 
presented.  For  convenience,  two  definitions  will  be  given  first. 

Definition  1 

Let  S  be  a  set  of  positive  integers  p^  (i  =  1,2 . s)  which  are 

smaller  than  or  equal  to  n,  the  order  of  the  matrices  A  and  B  of  the 
system  Ax  =  aBx,  i.e.,  S  =  (pj,  Pg.-.-.p  ).  Then,  R  is  defined  as  the 
subspace  spanned  by  the  s  eigenvectors  x^  (ieS). 

Definition  2 

If  no  vector  in  the  subspace  R  is  orthogonal  to  all  vectors 
Yj  (i  =  l,2,...,s)  (with  respect  to  B),  then  the  set  of  the  vectors 
yj  (i  =  l,2,...,s)  is  said  to  be  an  admissible  frame  with  respect  to  the 
subspace  R. 

Theorem  1 

With  the  above  definitions,  if  none  of  the  eigenvalues  A^  (ieS)  is 
equal  to  any  x.  (j^S),  then  among  all  admissible  frames  of  vectors 

J 

y-  (i  =  l,2,...,s)  a  frame  which  renders  w  extremum  in  the  following  con¬ 
strained  stationary- value  problem  spans  the  subspace  R,  and  its  stationary 
value  is  the  sum  of  the  eigenvalues  \ ^  (ieS): 
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Find  the  stationary  value  of 

s 

w  =  £  ylAy.  (C.l) 

i=l 

subject  to 

y^By..  =  6ij  (i,j  =  1,2,. . .  ,s)  (C.2) 

where  6. .  is  the  Kronecker  delta. 

1  J 

Proof:  For  convenience,  but  without  loss  of  generality,  the  set 
S  =  (l,2,...,s  n)  is  taken,  that  is,  R  is  the  subspace  spanned  by 
x.  (i  =  1,2, Let  the  vectors  y^  {i  =  1,2,. ...s)  be  expanded  in  a 
series  of  the  eigenvectors  xk  (k  *  l,2,...,n): 

n 

yi  =  £  Cki^k  (k  =  1,2,. . . ,n)  (C. 3) 

k-1 

It  will  be  shown  that  a  solution  of  Eqs.  (C.l)  and  (C.2)  yields 

ckl-  =  0  for  k  =  s+1,  s+2,...,n  (C. 5) 

Substitution  of  Eq.  (C.3)  into  Eqs.  (C.l)  and  (C.2)  results  in 
s  n 

W  =  £  £  xk  c*.  (C.6) 

i=l  i=l 


and 
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C._,  Cia  •  =  6,, 


k=l 


'ki  vkj  ij 


(i  »J  -  1,2,. . . ,s) 


(C.7) 


Let  us  use  the  method  of  Lagrange  multipliers  to  solve  the  stationary-value 
problem  of  Eqs.  (C.6)  and  (C.7).  Introducing  the  undetermined  multipliers 
v.-j  (i»j  =  1,2 . s)  and  letting  u . .  =  y  .  .,  we  have  the  Lagrangian 

•  J  I J  J  I 


s  n  s  s  n 

L  =  ^  ^  Akcki  L  uij  ckiCkj  '  <Sij) 

i=l  k=l  i=l  j=l  k=l 


(C.8) 


Since  the  first  derivatives  of  L  with  respect  to  the  unknowns  cki.  and 
should  vanish. 


9L 


3C 


ki 


=  2  (Akck.  -  =  0 

3=1 


(C.9) 


3L  _  r-  „  „  _  n 

kickj  "  3ij  "  u 


U 


k=l 


(C.10) 


We  may  write  Eqs.  (C.9)  and  (C.10)  in  matrix  form  as 


AC  =  CD 


CTC  =  I. 


where 


(C. 11) 

(C. 12) 


A  s  diag  (Xj»  ^2** '  *  *^n) 


I$  ■  unit  matrix  of  order  s 


and  the  (k,1)th  element  of  the  nxs  matrix  C  is  (k  =  l,2,...,n; 
i  =  1,2,...  ,s),  and  the  (i,j)th  element  of  the  sxs  matrix  D  is 
Mij  ^ Let  the  matrix  C  be  partitioned  into  the  two 
submatrices  Cg  and  Cn_$,  where  Cg  is  the  sxs  matrix  having  the  elements 
of  the  first  s  rows  of  C,  and  C  is  the  remaining  (n-s)xs  matrix. 
Then,  Eqs.  (C.ll)  and  (C.12)  may  be  written  as 


*,  Cs 

•CSD 

(C.13) 

Vs 

Cn-s 

■  C„-sD 

(C.14) 

ft  +  Cn-s 

Cn-s 

■  's 

(C. 15) 

where 

As 

®  diag  .. 

An-s 

=  diag  (xs+1*^s+2»- • • »xn) 

Since  y..  =  y..  was  taken 

’  J  J  * 

0T 

=  D 

(C. 16) 

by  which  from  Eq.  (C.13) 

CsAs 

=  DC] 

(C.17) 

Postmultiplication  of  Eq.  (C.14)  by  C  and  use  of  Eq.  (C.17)  leads  to 
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Let 


Then,  Eq.  (C.18)  yields 


(C. 19) 


Since 


Vsuij  =  Yu 


( i  =  1 ,2 , • . • ,n-s  i 
j  =  1,2,...  ,s)  (C.20) 


U  = 


(C.21) 


But  the  set  of  vectors  y  (i  =  1,2, ...,s)  is  an  admissible  frame  with 
respect  to  B,  from  which  it  is  not  difficult  to  show  that 


det  CJS  ?  0  (C.  22) 

From  Eqs.  (C.21)  and  (C.22),  we  obtain 

Vs  * 0  <c-23> 

or 

cki  —  0  (i  —  1,2 . s;  k  =  s+l,s+2,. . .  ,n)  (C.24) 

This  shows  that  the  subspace  spanned  by  the  vectors  (i  =  1,2 . s)  is 

the  subspace  of  the  eigenvectors  x.  (i  =  l,2,...,s),  which  is  to  be  proved 
here.  Furthermore,  from  Eq.  (C.6)  we  obtain 
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w 


Vki 


1=1  k=l 


s 


L. 


k=l 


s 


ki  =  l 


s 

-  Z  \  <c-25) 

k=l 

which  implies  that  the  stationary  value  is  the  sum  of  the  eigenvalues 
( k  =  1 *2 , . . . »s) . 

Theorem  2 

If  a  frame  of  s  vectors  y..  (i  =  l,2,...,s)  which  are  mutually  ortho¬ 
normal  with  respect  to  B  spans  the  subspace  R,  then  Eqs.  (C.13),  (C.14), 
and  (C.15)  are  satisfied,  i.e.,  w  is  stationary. 

Proof:  Since  the  vectors  y^  (i  =  l,2,...,s)  are  in  the  subspace  R, 
and  are  orthonormal  with  respect  to  B, 

Cn-s  B  0 
and 

clc  =  I  or  cl  =  C"1  (C.26) 

Hence,  Eqs.  (C.14)  and  (C.15)  are  satisfied.  Furthermore,  we  have 


(C.27) 


We  now  define  D  by  D  =  A$CS,  then 

AsCs  =  CsD  (C. 28) 

which  is  equivalent  to  Eq.  (C.13),  i.e.,  Eq.  (C.13)  is  also  satisfied. 
Theorem  3 

If  the  s  vectors  y^  (i  *  l,2,...,s)  span  the  subspace  R  of  x..  (ieS), 
then  the  Lagrange  multipliers  y. .  (i,j  =  l,2,...,s)  defined  in  Theorem  1 

'  J 

have  the  following  properties:  if  the  eigenvalues  x.  (ieS)  are  close 
together 

y^-  «  u^  for  i  f  j  (C.29) 

and  if  the  eigenvalues  are  multiple,  i.e..  A*  =  A^  (ieS) 

uid  =  C  for  i  t  j 

=  A*  (C.30) 

Proof:  For  convenience,  we  take  the  set  S  =  (l,2,...,s). 
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Then,  from  Eq.  (C.13) 


or 


From  Eq.  (C.7) 


D  -  Cs  Vs 


yij  =  £  xkckickj 
k=l 


(C. 31) 


=  x£  ckickj +  I  t^k  ■  V  ckickj  (c-32^ 


ik=l  k=l 


^  ckickj  =  6ij 
k=l 


(i » j  =  1,2,. . .  ,s) 


Thus 


"1J  *  £  <Ak  ‘  V  ckfckj  for  W  j 
k=l 


ii 


N  *  £  '  V  cm 

k=l 


(C. 33) 


If  the  eigenvalues  a ^  (i  =  1,2 . s)  are  close  together,  i.e. , 

|  Xk  "  Xi  |  «  X  ^  (k  1 ) ,  then  Eq.  (C.33)  implies  that 

|  yij  |  <<  |  yii  |  for  j  *  J  (C. 34) 


A-tm 
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Furthermore,  if  all  the  eigenvalues  x^  (i  =  l,2,...,s)  are  multiple,  i.e. 
x*  =  Xi  =  x2  = .  =  x  ,  then  from  Eq.  (C.33) 


p. .  =  0 


•Mi  =  xi  -  x* 


(i  =  1,2,. . . »s ) 


(C.35) 
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TABLE  1.  NUMBER  OF  OPERATIONS  FOR  EIGENSOLUTIONS 


Method 


Proposed 
Method  of 
Chapter  2 


Operation 


Calculation 


Multiplication  A  -  Aj^B 


Factorization 


LDU  =■  A  -  A.*0)B 


Number  of  Operations 


n  (mb  +  1) 

T  nma  (ma  +  3) 


Iteration 

Multiplication  Axj 


BXJ 


(k) 


(k) 


r  (k)  =  Ax/k)  -  AA.(k)Bx/k) 

J  J  J  J 


Factorization 


LDU  =  F 


00 


Solve  Eq.  (2.7)  for  Ax/k^  and  aa/1^ 


n  (2m  +  1) 

a 


n  (2mj}  +  1) 


n  (mfl  +  1) 


2n  (m  +  1) 

a 


Total 


N  =  i  pn  (in'"  +  3m  +  2m.  +  2)  +  T  n  (5m  +  2m.  +  6) 


Where 


-00  . 


A-A^0)B 


-x  (k)  B 

L  XJ  8 


-Bx 


00 


Tp  ■  Total  number  of  iterations  by  the  proposed  method-  Tp  '  Tr 

N  =  Total  number  of  operations  by  the  proposed  method. 
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TABLE  1.  (Continued) 


Method 

Operation 

Calculation 

Number  of  Operations 

Proposed 

Multiplication 

« v°'- 

LDU  =  A  - 

n  (m^  +  1) 

Method  of 
Chapter  3 

Factorization 

J  nma  (ma  +  3) 

Iteration 

Multiplication 

j 

n  (2ma  +  1) 

Mul tipi ication 

By/1'5  (i  =  1.2 . s) 

sn  (2^  +  1) 

Multiplicati  on 

-  *,“>  -  Z  .,j|k|  »y,(l) 

sn 

Factorization 

i=l 

LDU  =  F^k) 

sn  [ma  +  j  (s  +  1)] 

Solve  Eq.  (3.21)  for  and  ad..(k) 

n  (2m  +  s  +  1) 

a 

Total  Np  =  I 

pn  (  m|+3ma+2mh+2)  +  Ton  [(s+4)  ma+2sm| 

,  ♦  |  (s3+7s+4)] 
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TABLE  1 .  (Continued) 


Method 

Operation 

Calculation 

Number  of  Operations 

Robinson- 

Multiplication 

A  -  Vj^B 

n  (ma  +  1) 

Harris 

Method 

Multiplication 

=  <A  '  xj^k) 

n  (m  +  1) 

a 

Multiplication 

Bx.<k> 

n  (mb  +  1) 

Factorization 

LDU  =  G^k) 

A  nm  (m  +  5)  +  n 

C  d  d 

Solve  Eq.  (2.6)  fo 

ir  ax.^  and  aa.^ 

J  J 

2n  (m  +  1) 

a 

Total 

Nr  =  ?  Trn  {ma  + 

13m  +  6  .  +  12) 
a  mb 

Where 


,00  - 


A-Vj^B 


T 

N 


r 

r 


=  Total  number  of  iterations  by  the  Robinson-Harris  method.  T  < 
-  Total  number  of  operations  by  the  Robinson-Harris  method. 


TABLE  1.  (Continued) 


Method 

Operation 

Calculation 

Number  of  Operations 

Subspace 

Factorization 

LOU  =■  A 

nm.  {m.  +  3)/2 

Iteration 

Method 

Iteration 

o  a 

Multiplication 

BX*k_1) 

qn  (2mb  +  1) 

Solve  for  Y*k) 

Ay<k)  =  bx^*1) 

qn  (2ma  +  1) 

Multiplication 

A*(k)  =  Y00TBX(k-i) 

qn  (q  +  l)/2 

Multiplication 

BY< k ) 

qn  (2mb  +  1) 

Multiplication 

B*(k)  M  y(k)TBY(k) 

qn  (q  +  l)/2 

Solve  for  Z(k) 

A*(k)z(k)  „  B*(k)z(k)D{k) 

0  (q^)  neglected 

and  D(k) 

Multiplication 

x(k)  .  y(k)z(k) 

nq2 

Subtotal 

qn  (2ma  +  411^  +  2q  +  4) 

Sturm  Sequence 

Check 

Multiplication 

A-\p(k)B 

n  (mfa  +  1) 

Factorization 

LOU  =  A-Ap^k)B 

nm,  (m  +  3)/2 

a  a 

Total  N 

'$  =  Tsqn  ^a*  +  4mb  +  2q  +  4) 

2 

+  n  (m  +  3m,  +  m.  +  1) 
a  a  o 

Note:  q 

-  max  (2p,  p+8) 

Ts 

=  Total  number  of 

iterations  by  subspace  iteration 

method. 

Ns 

-  Total  number  of  ooerations  by  subsoace  iteration 

method. 

It 

is  assumed  that  m 

a  -  V 
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TABLE  2.  EIGENVALUES  OF  THE  PLANE  FRAME  PROBLEM 
(DISTINCT  ROOTS) 


Method 

of 

Analysis 

Iteration 

Number 

(1) 

Eigenvalues 

(2)  (3) 

(4) 

Proposed 

0 

0.474744x10° 

0.443880X101 

0.132953xl02 

0. 284745xl02 

Method 

(0.36xl0-2)* 

(0.28X10*1) 

(0.77x10**) 

(0.22x10°) 

1 

0.474744x10° 

0.443876X101 

0.132921xl02 

0.284091xl02 

(0.44xl0*5) 

(0.45xl0*3) 

(0.38xl0*2) 

(0.27x10**) 

2 

0.474744x10° 

0.443876x10* 

0.132921xl02 

0.284091xl02 

(0.82xl0*13) 

(0.27x10*®) 

(0. 17xl0*6) 

(0. 18xl0*3) 

3 

0.284091xl02 
(0. 19x10*®) 

4 

0.284091xl02 

(0.23X10*7) 

Robinson 

0 

0.474744x10° 

0.443880x10* 

0. 132953xl02 

0.284745xl02 

Harris 

Method 

(0.36xl0*2) 

(0.28x10**) 

(0.77x10**) 

(0.22x10° 

1 

0.474744x10° 

0.443876x10* 

0. 132921xl02 

0.284091xl02 

(0.44x10*®) 

(0. 45xl0*3) 

(0. 38xl0*2) 

(0.27x10**) 

2 

0.474744x10° 

(0.82X10*13) 

0.443876x10* 

(0.27xl0*9) 

0. 132921xl02 
(0. 17x10*®) 

0.284091xl02 
(0. 18xl0*3) 

3 

0.284091xl02 
(0. 13xl0*8) 

*  :  Numbers  In  parentheses  Indicate  errors  In  the  approximate 
(k) 

eigenvectors  xy 


* 
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TABLE  2.  (Continued) 


Method 

of 

Analysis 

Iteration 

Number 

(1) 

Eigenvalues 

(2)  (3) 

(4) 

Subspace 

Iteration 

I 

0.476915x10° 

n 

0.465927X101 

0.443880X101 

(0.28X10'1) 

0.153636xl02 

0.369761xl02 

Method 

2 

0.474744x10° 

<0.36xl0”2)* 

0. 132953x10* 
(0. 77x10” 1 ) 

0.284745x10* 

(0.22xl02) 

3 

0.474744x10° 

(0.44x10*°) 

0.443876X101 

(0.44xl0'3) 

0. 132921xl02 
(0.35xl0‘2) 

0.284116xl02 

(0.19X10'1) 

4 

0.474744x10° 
(0. llxlO"7) 

0.443876X101 

(0.22xl0”4) 

0. 132921xl02 
(0.33xl0”3) 

0. 284099xl02 
(0.48xl0‘2) 

5 

0.474744x10° 

(0.83xl0”10) 

0.443876X101 
(0. 19xl0"5) 

0. 13292 lxlO2 
(0. 90xl0'4) 

0.284094xl02 
(0. 34xl02) 

6 

0.474744x10° 

(0.71xl0”12) 

0.443876x10* 
(0. 16x10"°) 

0.132921xl02 

(0.28xl0~4) 

0.284092xl02 

(0.36xl02) 

7 

0.474744x10° 

(O.llxlO*13) 

0.443876X101 
(0. 57x10”°) 

0. 132921xl02 
(0. 28x10”°) 

0.284091X102 
(0. llxlO”2) 

8 

0.474744x10° 

(0.13X10*13) 

0.443876X101 
(0.  Z6xl0-9) 

0. 132921xl02 
(0. 28x10°) 

0.284091xl02 
(0. 14xl03) 

9 

0.474744x10° 
(0. 16xl0"13) 

0.443876X101 
(0. 17x10” 10) 

0.132921xl02 

(0.53xl0'7) 

0.284091xl02 

(0.32xl0’4) 

10 

0.474744x10° 

(0.13xl0*13) 

0.443876X101 

(0.93xl0”12) 

0.1 3292 lx 102 
(0. 91x10°) 

0.284091xl02 

(0.99x10°) 

11 

0.474744x10° 

0.443876x10* 

0. 13292 lxlO2 

0.284091xl02 

(0. 13x10' 13) 

(0.52xl0”13) 

(0.12x10”°) 

(0.33x10”°) 

12 

0.474744x10° 

(0.28X10'13) 

0.443876x10* 

(0.50xl0'13) 

0. 132921xl02 
(0. 18x10°) 

0.284091xl02 

(0.12x10”°) 

*  :  Numbers  in  parentheses  Indicate  errors  in  the  approximate 
eigenvectors  x/ 


TABLE  3.  EIGENVALUES  OF  THE  CIRCULAR  ARCH  PROBLEM 
(DISTINCT  ROOTS) 


Method 

of 

Analysis 

Iteration 

Number 

(1) 

Eigenvalues/)^* 

(2) 

(3) 

Proposed 

0 

0.102714xl0'2 

0.919516xl0”2 

0.363073X10*1 

Method 

1 

0.102640xl0'2 

0.909467xl0”2 

0.341422X10”1 

2 

(0.59xl0"3)** 

0. 102640xl0*2 

(0. 53xl0”2) 
0.909467xl0”2 

(0.35X10'1) 

0.341448X10”1 

3 

4 

(0.26x10”®) 

(O.lOxlO'5) 

0.909467xl0*2 

(O.lOxlO”8) 

(0.21xl0‘3) 

0.34144SX10”1 

(0.17xl0'5) 

0.341448X10'1 
(0. 33x10* 7 ) 

Robinson-  0 

Harris 

0.102714xl0”2 

0.919516xl0”2 

O 

0.363073x10 

Method  1 

0.102640x10  c 

0.909467x10'^ 

0.341422x10' 

(0.59xl0*3) 

(0.53xl0'2) 

(0.35X10'1) 

2 

0. 102640xl0'2 

0.909467xl0"2 

0.341443x10' 

(0.26x10'®) 

(O.lOxlO'5) 

(0. 21xl0'3) 

3 

0.909467xl0'2 

(0.45xl0'12) 

0.341443x10' 

(0.20x10'®) 

*  :  *0  *  E/Pa2(l- 

v2). 

:  Numbers  In  parentheses  Indicate  errors  in  the  approximate 
eigenvectors  x^  . 


y 

b 


TABLE  3.  (Continued) 


j 

4  -« 

Method 

Eigenvalues/ \n* 

of 

Iteration 

Analysis 

Number 

(1) 

(2) 

(3) 

«  * 

Subspace 

1 

0. 102714xl0*2 

0.919516X10"2 

0.363073X10"1 

Iteration 

0.341476X10"1 

Method 

2 

0.102640x10  c 

0.909468x10** 

(0.59xl0"3)** 

(0.52xl0‘2) 

(0.32X10"1) 

3 

0. 102640xl0"2 

0.909467xl0*2 

0. 341448x10* 1 

(0.41x10*®) 

(0.62x10*®) 

(0.26xl0"2) 

•  • 

4 

0. 102640xl0"2 

0.909467xl0"2 

0.341448X10*1 

(0.58x10"®) 

(0. 12xl0*5) 

(0.23xl0*3) 

5 

0. 102640xl0”2 

0. 909467xl0*Z 

0.341448X10*1 

* 

(0. 12xl0~ 11 ) 

(0.28xl0‘7) 

(0.21x10"®) 

•  - 

6 

0. 102640xl0*2 

0.909467xl0"2 

0.341448X10*1 

«  • 

(0.23xl0*13) 

(0.79x10*®) 

(0.19xl0‘5) 

7 

0. 102640xl0*2 

0.909467xl0"2 

0. 341448x10" 1 

(0. 12x10" l3) 

(0.25x10" 30) 

(O.PxlO*6) 

•  — 

j 

*  * 


*  :  *Q  *  E/Ba2(l-v2). 

**  :  Numbers  in  parentheses  indicate  errors  in  the  approximate 
(k) 

eigenvectors  x.'  . 

I 

l 

I 

I 

I 

I 
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TABLE  4.  EIGENVALUES  OF  THE  SQUARE  PLATE  PROBLEM 
(DOUBLE  ROOTS) 


Method 

of 

Analysis 

Iteration 

Number 

(1) 

Eigenvalues/a* 

(2)  (3) 

(4) 

Proposed 

0 

0.375840x10* 

0.2302 12xl02 

0.230212xl02 

0. 530627xl02 

Method 

I 

0.375838x10* 

(0.80xl0’3)** 

0.230121xl02 

(0.10x10’*) 

0. 23012  lx 102 
(0.10x10’*) 

0.529932xl02 

(0.21x10**) 

2 

0.375838x10* 

(0.44x10’**) 

0.230121xl02 

(0.46xl0*7) 

0.230121xl02 

(0.46xl0’7) 

0. 529932xl02 
(0. 59xl0’6) 

Subspace 

Iteration 

1 

0.375840x10* 

0.375838x10* 

(0.80xl0*3) 

0.230212X102 

0.230212xl02 

0. 53062 7xl02 

Method 

2 

0.230121x10* 

(0.10x10**) 

0.230121x10“' 

(0.10x10"*) 

0. 529932xl0<: 
(0.21x10**) 

3 

0.375838x10* 

(0.65X10"6) 

0.230121xl02 

(O.llxlO’3) 

0.230121xl02 

(O.llxlO'3) 

0. 529932xl02 
(0. 43xl0'3) 

4 

0.375838x10* 
(0.53x 10"9) 

0.230121x10^ 

(0.13xl0"5) 

0.230121xl02 

(0.13X10*5) 

0. 529932xl02 
(0.90xl0*5) 

5 

0.375838x10* 
(0. 44x10* *2) 

0.230121xl02 
(0. 1 5x 10’ 7 ) 

0.230121xl02 
(0. 15xlO’7) 

0. 529932xl02 
(0. 19xl0*6) 

*  :  a  =  ir4De/(a4p).  where  Dg  *  Eh3/12(l  -  u2). 

**  :  Numbers  in  parentheses  indicate  errors  in  the  approximate 
(k) 

eigenvectors  yj 


TABLE  5.  EIGENVALUES  OF  THE  RECTANGULAR  PLATE  PROBLEM 

(CLOSE  ROOTS) 


Method 

of 

Analysis 

Iteration 

Number 

(1) 

El gen values/a* 

(2)  (3) 

(4) 

Proposed 

0 

0.368470X101 

0. 222957xl02 

0. 228454xl02 

0.520215xl02 

Method 

1 

0.268468X101 

(0.80xl0*3) 

0.222868xl02 

(O.lOxlO'1) 

0.228264xl02 

(O.lOxlO*1) 

0. 519522xl02 
(0. 21x10**) 

2 

0. 368468x10* 
(0.44x10' 11 ) 

0. 222868xl02 
(0. 47xl0'7) 

0.228364xl02 

(0.45xl0*7) 

0.519533xl02 

(0.59xl0*6) 

Subspace 

Iteration 

1 

0. 368470x10* 

0. 368468x10* 
(0.80xl0'3) 

0.222957xl02 

0.228454xl02 

0. 520215xl02 

Method 

2 

0.222868x10*" 

(O.lOxlO'1) 

0.228364x10* 

(O.lOxlO*1) 

0.519533x10* 

(0.21x10**) 

3 

0.368468X101 

(0.65xl0'6) 

0.222868xl02 

(O.llxlO'3) 

0.228364xl02 

(O.llxlO*3) 

0. 519533xl02 
(0.43xl0*3) 

4 

0.368468X101 

(0.53x10"®) 

0.222868xl02 
(0. 13x10*®) 

0. 228364xl0Z 
(0. 13x10*®) 

0.519533xl02 

(0.90x10*®) 

5 

0. 368468x10* 
(0.46xl0'12) 

0.222868xl02 

(0.15xl0'7) 

0.228364xl02 
(0. 14xlO*7) 

0. 519533xl02 
(0. 19x10*®) 

*  :  a  =  *4De/(a4p),  where  Dg  »  Eh3/12(l  -  v2). 

**  :  Numbers  in  parentheses  indicate  errors  in  the  approximate 
(k) 

eigenvectors  y^ 
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TABLE  6.  COMPARISON  OF  THE  TOTAL  NUMBER  OF  OPERATIONS 


Input  Data  Number  of  Number  of  Operations  Ratio 

Problem  Iterations 


Type 

n  ma  nip  p  q 

TP 

Tr 

Ts 

"p 

Nr 

Ns 

VNp  VN, 

Frame 

330  35  35  4  10 

10 

9 

12 

3.50x10® 

4.57x10® 

9.72x10® 

1.31  2.78 

Arch 

22  4  0  3  5 

9 

8 

7 

8.87xl03 

9.77xl03 

1.76xl04 

1.10  1.98 

Plate 

39  16  16  4  9 

8 

- 

5 

1.27x10® 

- 

2.20x10® 

1.73 

Note 

n  :  Order  of  stiffness  and  mass  matrices 

in  :  Average  half  bandwidth  of  stiffness  matrix 

d 

mb  :  Average  half  bandwidth  of  mass  matrix 

p  :  Number  of  eigenvalues  and  eigenvectors  sought 

q  :  Number  of  iteration  vectors,  q  =  max  (2p,  p+8) 

Tp  :  Number  of  iterations  by  the  proposed  method 

Tr  :  Number  of  iterations  by  the  Robinson-Harris  method 

T$  :  Number  of  iterations  by  the  subspace  iteration  method 

Np  :  Total  number  of  operations  by  the  proposed  method 

Nr  :  Total  number  of  operations  by  the  Robinson-Harris  method 

Ns  :  Total  number  of  operations  by  the  subspace  iteration  method 
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TABLE  7.  COMPARISON  BETWEEN  THE  THEORETICAL  CONVERGENCE  RATES 
FOR  EIGENVECTORS  AND  THE  NUMERICAL  RESULTS  -  FRAME  PROBLEM 

(DISTINCT  ROOTS) 


Method 

Eigenvalue  Number 

of 

Iteration 

Analysis 

Number 

1 

2 

3 

4 

Proposed 

1 

1.9x10*® 

6. 0x10* 7 

4.4x10*® 

6. 7xl03 

Method 

2 

1.1x10'® 

3 

1.3x10"® 

1.2x10'® 

Theory 

1.0x10*® 

3.6xl0'4 

1.2x10'® 

Subspace 

2 

1.2xl0*3 

1.6x10*® 

4.5x10'® 

8.6x10'® 

Iteration 

Method 

3 

4 

2.5xl0*3 

7.5xl0*3 

5.0x10*® 

8.7x10'® 

9.4x10'® 

2.7X10"1 

2.5xl0'3 
7. lxlO*1 

5 

8.6xl0"3 

8.4x10*® 

3. 1x10' 1 

1.1x10"° 

6 

1. 5xl0*3 

3.6x10*® 

l.OxlO"1 

3. 1x10' 1 

7 

* 

4.6x10*® 

l.OxlO'1 

1.3xl03 

8 

* 

6.6x10*® 

1.9X10'1 

2. 3x10' 1 

9 

* 

5.5x10"® 

1.7X10'1 

3. lxlO3 

10 

* 

5.6x10*® 

1.3x10" 3 

3.3xl0'3 

11 

* 

* 

l.SxlO'1 

3.6xl0'3 

Theory 

5.6xl0*3 

5.2x10*® 

1.6x10* 3 

3.3xlO'3 

*  :  Errors  too  small  for  comparison  because  of  round-off  error. 
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TABLE  8.  COMPARISON  BETWEEN  THE  THEORETICAL  CONVERGENCE  RATES  FOR 
EIGENVECTORS  AND  THE  NUMERICAL  RESULTS  -  SQUARE  PLATE  PROBLEM 

(DOUBLE  ROOTS) 


Method 

Eigenvalue  Number 

of 

Iteration 

Analysts 

Number 

1 

2 

3 

4 

Proposed 

1 

5.5x10'® 

4.6x10'® 

4.6x10'® 

2.8x10'® 

Method 

Theory 

l.OxlO'6 

4.7xl0-4 

4.7xl0'4 

2. 3x10' 3 

Subspace 

2 

8.1xl0'4 

l.lxlO'2 

l.lxlO'2 

2.0xl0*2 

Iteration 

Method 

3 

4 

8.2xl0'4 

8.3xl0~4 

1.2x10’ 2 
1.2xl0'2 

1.2xl0'2 

1.2xl0-2 

2.1xl0‘2 
2.  lxlO"2 

Theory 

l.lxlO'2 

6.8xl0'2 

6.8xl0'2 

1.6X10*1 
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TABLE  9.  NUMERICAL  CONVERGENCE  RATES  FOR  EIGENVECTORS  - 
RECTANGULAR  PLATE  PROBLEM 
(CLOSE  ROOTS) 


Method 

of 

Analysis 

Iteration 

Number 

1 

Eigenvalue  Number 

2  3 

4 

Proposed 

Method 

1 

5.5xl0'9 

4.7x10"® 

4. 5xl0"6 

2.8xl0"5 

Subs pc  ;e 

2 

8. lxlO'4 

l.lxlO"2 

l.lxlO'2 

2.0xl02 

Iteration 

Method 

3 

8. 2xl0'4 

1.2xl0"2 

1.2xlO*2 

2.  lxlO2 

4 

8. 7xlO"4 

1.2xl0"2 

l.lxlO"2 

2. lxlO"2 

Theory 

l.lxlO'2 

6. 8xl0’2 

7. OxlO*2 

1.6X10"1 

r 


For  All  Beams  and  Columns 

Area  of  Cross-Section  A 
Moment  of  Inertia  of  Cross-Section  I 
Young's  Modulus  E 
Mass  Density  p 


=  2. 787x1 0-1  m2 
=  8.631xl0"3  m4 
=  2. 068x1 010  Pa 
=  1.602xl04  kg/m3 


FIG.  2  TEN-STORY,  TEN-BAY  PLANE  FRAME 
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Navy  (Con't . ) 

Naval  Research  Laboratory 
Washington,  D.C.  20375 
Attn:  Code  8400 
8410 
8430 
8440 
6300 
6390 
6380 

David  W.  Taylor  Naval  Ship  Research 
and  Development  Center 
Annapolis,  Maryland  21402 
Attn:  Code  2740 
28 
281 

Naval  Weapons  Center 
China  Lake,  California  93555 
Attn:  Code  4062 
4520 

Comsianding  Officer 

Naval  Civil  Engineering  Laboratory 

Code  L31 

Port  Hueneme,  California  93041 

Naval  Surface  Weapons  Center 
White  Oak 

Silver  Spring,  Maryland  20910 
Attn:  Code  R-10 
G-402 
K-82 

Technical  Director 

Naval  Ocean  Systems  Center 

San  Diego,  California  92152 

Supervisor  of  Shipbuilding 
U.S.  Navy 

Newport  News,  Virginia  23607 

Navy  Underwater  Sound 
Reference  Division 
Naval  Research  Laboratory 
P.0.  Box  8337 
Orlando,  Florida  32806 

Chief  of  Naval  Operations 
Department  of  the  Navy 
Washington,  D.C.  20350 
Attn:  Code  OP-098 
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Navy  (Con't.) 


Strategic  Systems  Project  Office 
Department  of  the  Navy 
Washington,  D.C.  20376 
Attn:  NS P-2 00 

Naval  Air  Systems  Command 
Department  of  the  Navy 
Washington,  D.C.  20361 

Attn:  Code  5302  (Aerospace  and  Structures) 
604  (Technical  Library) 

320B  (Structures) 

Naval  Air  Development  Center 
Warminster,  Pennsylvania  18974 
Attn:  Aerospace  Mechanics 
Code  606 

U.S.  Naval  Academy 
Engineering  Department 
Annapolis,  Maryland  21402 

Naval  Facilities  Engineering  Command 
200  Stovall  Street 
Alexandria,  Virginia  22332 
Attn:  Code  03  (Research  and  Development) 
04B 
045 

14114  (Technical  Library) 

Naval  Sea  Systems  Command 
Department  of  the  Navy 
Washington,  D.C.  20362 
Attn:  Code  05H 
312 

322 

323 
05R 
32R 


Commander  and  Director 
David  W.  Taylor  Naval  Ship 

Research  and  Development  Center 
Bethesda,  Maryland  20084 
Attn:  Code  042 
17 

172 

173 

174 
1800 
1844 

012.2 

1900 

1901 
1945 
1960 
1962 

Naval  Underwater  Systems  Center 
Newport ,  Rhode  Island  02840 
Attn:  Bruce  Sandman,  Code  3634 

Naval  Surface  Weapons  Center 
Dahlgren  Laboratory 
Dahlgren,  Virginia  22448 
Attn:  Code  G04 
G20 

Technical  Director 

Mare  Island  Naval  Shipyard 

Vallejo,  California  94592 

U.S.  Naval  Postgraduate  School 

Library 

Code  0384 

Monterey,  California  93940 

Webb  Institute  of  Naval  Architecture 
Attn:  Librarian 
Crescent  Beach  Road,  Glen  Cove 
Long  Island,  New  York  11542 


Commanding  Officer  (2) 

U.S.  Army  Research  Office 
P.0.  Box  12211 

Research  Triangle  Park,  NC  27709 
Attn:  Mr.  J.  J.  Murray,  CRD-AA-IP 
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Army  (Con' t. ) 

Watervliet  Arsenal 
MAGGS  Research  Center 
Watervliet,  New  York  12189 
Attn:  Director  of  Research 

U.S.  Army  Materials  and  Mechanics 
Research  Center 

Watertown,  Massachusetts  02172 
Attn:  Dr.  R.  Shea,  DRXMR-T 

U.S.  Army  Missile  Research  and 
Development  Center 
Redstone  Scientific  Information 
Center 

Chief,  Document  Section 
Redstone  Arsenal,  Alabama  3S809 

Army  Research  and  Development 
Center 

Fort  Belvoir,  Virginia  22060 
NASA 

National  Aeronautics  and  Space 
Administration 
Structures  Research  Division 
Langley  Research  Center 
Langley  Station 
Hampton,  Virginia  23365 

National  Aeronautics  and  Space 
Adminia  tration 

Associate  Administrator  for  Advanced 
Research  and  Technology 
Washington,  D.C.  20546 

Air  Force 

Wright-Patterson  Air  Force  Base 
Dayton,  Ohio  45433 
Attn:  AFFDL  (FB) 

(FBR) 

(FBE) 

(FBS) 

AFML  (MBM) 

Chief  Applied  Mechanics  Group 
U.S.  Air  Force  Institute  of  Technology 
Wright-Patterson  Air  Force  Base 
Dayton,  Ohio  45433 


Air  Force  (Con't.) 

Chief,  Civil  Engineering  Branch 
WLRC,  Research  Division 
Air  Force  Weapons  Laboratory 
Kirtland  Air  Force  Base 
Albuquerque,  New  Mexico  87117 

Air  Force  Office  of  Scientific  Research 
Bolling  Air  Force  Base 
Washington,  D.C.  20332 
Attn:  Mechanics  Division 

Department  of  the  Air  Force 
Air  University  Library 
Maxwell  Air  Force  Base 
Montgomery,  Alabama  36112 

Other  Government  Activities 

Commandant 

Chief,  Testing  and  Development  Division 
U.S.  Coast  Guard 
1300  E  Street,  NW. 

Washington,  D.C.  20226 

Technical  Director 
Marine  Corps  Development 
and  Education  Command 
Quantico,  Virginia  22134 

Director  Defense  Research 
and  Engineering 
Technical  Library 
Room  3C128 
The  Pentagon 
Washington,  D.C.  20301 

Dr.  M.  Gaus 

National  Science  Foundation 
Environmental  Research  Division 
Washington,  D.C.  20550 

Library  of  Congress 

Science  and  Technology  Division 

Washington,  D.C.  20540 

Director 

Defense  Nuclear  Agency 
Washington,  D.C.  20305 
Attn:  SPSS 
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Other  Government  Activities  (Con't) 

Mr.  Jerome  Persh 
Staff  Specialist  for  Materials 
and  Structures 
OUSDR&E ,  The  Pentagon 
Room  3D1089 

Washington,  D.C.  20301 

Chief,  Airframe  and  Equipment  Branch 
FS-120 

Office  of  Flight  Standards 
Federal  Aviation  Agency 
Washington,  D.C.  20553 

National  Academy  of  Sciences 
National  Research  Council 
Ship  Hull  Research  Committee 
2101  Constitution  Avenue 
Washington,  D.C.  20418 
Attn:  Mr.  A.  R.  Lytle 

National  Science  Foundation 
Engineering  Mechanics  Section 
Division  of  Engineering 
Washington,  D.C.  20550 

Picatinny  Arsenal 

Plastics  Technical  Evaluation  Center 
Attn:  Technical  Information  Section 
Dover,  New  Jersey  07801 

Maritime  Administration 
Office  of  Maritime  Technology 
14th  and  Constitution  Avenue,  NW. 
Washington,  D.C.  20230 

PART  2  -  Contractors  and  Other  Technical 
Collaborators 

Universities 

Dr.  J.  Tinsley  Oden 
University  of  Texas  at  Austin 
345  Engineering  Science  Building 
Austin,  Texas  78712 

Professor  Julius  Miklowitz 
California  Institute  of  Technology 
Division  of  Engineering 
and  Applied  Sciences 
Pasadena,  California  91109 


Universities  (Con't) 

Dr.  Harold  Liebowitz,  Dean 
school  of  Engineering  and 
Applied  Science 
George  Washington  University 
Washington,  D.C.  20052 

Professor  Eli  Sternberg 
California  Institute  of  Technology 
Division  of  Engineering  and 
Applied  Sciences 
Pasadena,  California  91109 

Professor  Paul  M.  Naghdi 
University  of  California 
Department  of  Mechanical  Engineering 
Berkeley,  California  94720 

Professor  A.  J.  Durelli 
Oakland  University 
School  of  Engineering 
Rochester,  Missouri  48063 

Professor  F.  L.  DiMaggio 
Columbia  University 
Department  of  Civil  Engineering 
New  York,  New  York  10027 

Professor  Norman  Jones 

The  University  of  Liverpool 

Department  of  Mechanical  Engineering 

P.  0.  Box  147 

Brownlow  Hill 

Liverpool  L69  3BX 

England 

Professor  E.  J.  Skudrzyk 
Pennsylvania  State  University 
Applied  Research  Laboratory 
Department  of  Physics 
State  College,  Pennsylvania  16801 

Professor  J.  Klosner 
Polytechnic  Institute  of  New  York 
Department  of  Mechanical  and 
Aerospace  Engineering 
333  Jay  Street 
Brooklyn,  New  York  11201 

Professor  R.  A.  Sch apery 
Texas  A&M  University 
Department  of  Civil  Engineering 
College  Station,  Texas  77843 
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Universities  (Con't.) 

Professor  Walter  D.  Pilkey 
University  of  Virginia 
Research  Laboratories  for  the 
Engineering  Sciences  and 
Applied  Sciences 
Charlottesville,  Virginia  22901 

Professor  K.  D.  WiL Inert 
Clarkson  College  of  Technology 
Department  of  Mechanical  Engineering 
Potsdam,  New  York  13676 

Dr.  Walter  E.  Haisler 
Texas  A&M  University 
Aerospace  Engineering  Department 
College  Station,  Texas  77843 

Dr.  Hussein  A.  Kamel 
University  of  Arizona 
Department  of  Aerospace  and 
Mechanical  Engineering 
Tucson,  Arizona  85721 

Dr.  S.  J.  Fenvea 
Carnegie-Mellon  University 
Department  of  Civil  Engineering 
Schenley  Park 

Pittsburgh,  Pennsylvania  15213 

Dr.  Ronald  L.  Huston 
Department  of  Engineering  Analysis 
University  of  Cincinnati 
Cincinnati,  Ohio  45221 

Professor  G.  C.  M.  Sih 
Lehigh  University 
Institute  of  Fracture  and 
Solid  Mechanics 
Bethlehem,  Pennsylvania  18U15 

Professor  Albert  S.  Kobayashi 
University  of  Washington 
Department  of  Mechanical  Engineering 
Seattle,  Washington  98105 

Professor  Daniel  Frederick 
Virginia  Polytechnic  Institute  and 
State  University 

Department  of  Engineering  Mechanics 
Blacksburg,  Virginia  24061 


Universities  (Con't) 

Professor  A.  C.  Eringen 
Princeton  University 
Department  of  Aerospace  and 
Mechanical  Sciences 
Princeton,  New  Jersey  08540 

Professor  E.  H.  Lee 

Stanford  University 

Division  of  Engineering  Mechanics 

Stanford,  California  94305 

Professor  Albert  I.  King 
Wayne  State  University 
Biomechanics  Research  Center 
Detroit,  Michigan  48202 

Dr.  V.  R.  Hodgson 
Wayne  State  University 
School  of  Medicine 
Detroit,  Michigan  48202 

Dean  B.  A.  Bo ley 
Northwestern  University 
Department  of  Civil  Engineering 
Evanston,  Illinois  60201 

Professor  P.  G.  Hodge,  Jr. 
University  of  Minnesota 
Department  of  Aerospace  Engineering 
and  Mechanics 

Minneapolis,  Minnesota  55455 

Dr.  D.  C.  Drucker 
University  of  Illinois 
Dean  of  Engineering 
Urbana,  Illinois  61801 

Professor  N.  M.  Newmark 
University  of  Illinois 
Department  of  Civil  Engineering 
Urbana,  Illinois  61803 

Professor  E.  Reissner 
University  of  California,  San  Diego 
Department  of  Applied  Mechanics 
La  Jolla,  California  92037 

Professor  William  A.  Nash 
University  of  Massachusetts 
Department  of  Mechanics  and 
Aerospace  Engineering 
Amherst,  Massachusetts  01002 
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Universit ies  (Con't) 

Professor  G.  Herrmann 
Stanford  University 
Department  of  Applied  Mechanics 
Stanford,  California  9430S 

Professor  J.  D.  Achenbach 
Northwest  University 
Department  of  Civil  Engineering 
Evanston,  Illinois  60201 

Professor  S.  B.  Dong 
University  of  California 
Department  of  Mechanics 
Los  Angeles,  California  90024 

Professor  Burt  Paul 
University  of  Pennsylvania 
Towne  School  of  Civil  and 
Mechanical  Engineering 
Philadelphia,  Pennsylvania  19104 

Professor  H.  W.  Liu 
Syracuse  University 
Department  of  Chemical  Engineering 
and  Metallurgy 
Syracuse,  New  York  13210 

Professor  S.  Bodner 
Technion  R&D  Foundation 
Haifa,  Israel 

Professor  Werner  Goldsmith 
University  of  California 
Department  of  Mechanical  Engineering 
Berkeley,  California  94720 

Professor  R.  S.  Rivlin 
Lehigh  University 
Center  for  the  Application 
of  Mathematics 

Bethlehem,  Pennsylvania  18015 

Professor  F.  A.  Cozzarelli 
State  University  of  New  York  at 
Buffalo 

Division  of  Interdisciplinary  Studies 
Karr  Parker  Engineering  Building 
Chemistry  Road 
Buffalo,  New  York  14214 


Universities  (Con't) 

Professor  Joseph  L.  Rose 
Drexel  University 

Department  of  Mechanical  Engineering 
and  Mechanics 

Philadelphia,  Pennsylvania  19104 

Professor  B.  K.  Donaldson 
University  of  Maryland 
Aerospace  Engineering  Department 
College  Park,  Maryland  20742 

Professor  Joseph  A.  Clark 
Catholic  University  of  America 
Department  of  Mechanical  Engineering 
Washington,  D.C.  20064 

Dr.  Samuel  B.  Batdorf 
University  of  California 
School  of  Engineering 
and  Applied  Science 
Los  Angeles,  California  90024 

Professor  Isaac  Fried 
Boston  University 
Department  of  Mathematics 
Boston,  Massachusetts  02215 

Professor  E.  Krempl 
Rensselaer  Polytechnic  Institute 
Division  of  Engineering 
Engineering  Mechanics 
Troy,  New  York  12181 

Dr.  Jack  R.  Vinson 
University  of  Delaware 
Department  of  Mechanical  and  Aerospace 
Engineering  and  the  Center  for 
Composite  Materials 
Newark,  Delaware  19711 

Dr .  J .  Duffy 
Brown  University 
Division  of  Engineering 
Providence,  Rhode  Island  02912 

Dr.  J.  L.  Swedlow 
Carnegie-Mellon  University 
Department  of  Mechanical  Engineering 
Pittsburgh,  Pennsylvania  15213 
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Universities  (Con't)  Universities  (Con't) 

Dr.  V.  K.  Varadan  Professor  V.  H.  Neubert 

Ohio  State  University  Research  Foundation  Pennsylvania  State  University 


Department  of  Engineering  Mechanics 
Columbus,  Ohio  43210 

Dr.  Z.  Hashin 

University  of  Pennsylvania 
Department  of  Metallurgy  and 
Materials  Science 
College  of  Engineering  and 
Applied  Science 

Philadelphia,  Pennsylvania  19104 

Dr.  Jackson  C.  S.  Yang 
University  of  Maryland 
Department  of  Mechanical  Engineering 
College  Park,  Maryland  20742 

Professor  T.  Y.  Chang 
University  of  Akron 
Department  of  Civil  Engineering 
Akron,  Ohio  44325 

Professor  Charles  W.  Bert 
University  of  Oklahoma 
School  of  Aerospace,  Mechanical, 
and  Nuclear  Engineering 
Norman,  Oklahoma  73019 

Professor  Satya  N.  Atluri 
Georgia  Institute  of  Technology 
School  of  Engineering  and 
Mechanics 

Atlanta,  Georgia  30332 

Professor  Graham  F.  Carey 
University  of  Texas  at  Austin 
Department  of  Aerospace  Engineering 
and  Engineering  Mechanics 
Austin,  Texas  78712 

Dr .  S .  S .  Wang 
University  of  Illinois 
Department  of  Theoretical  and 
Applied  Mechanics 
Urbana,  Illinois  61801 

Professor  J.  F.  Abel 
Cornell  University 
Department  of  Theoretical 
and  Applied  Mechanics 
Ithaca,  New  York  14853 


Department  of  Engineering  Science 
and  Mechanics 

University  Park,  Pennsylvania  16802 

Professor  A.  W.  Leissa 
Ohio  State  University 
Department  of  Engineering  Mechanics 
Columbus,  Ohio  43212 

Professor  C.  A.  Brebbia 
University  of  California,  Irvine 
Department  of  Civil  Engineering 
School  of  Engineering 
Irvine,  California  92717 

Dr.  George  T.  Hahn 
Vanderbilt  University 
Mechanical  Engineering  and 
Materials  Science 
Nashville,  Tennessee  37235 

Dean  Richard  H.  Gallagher 
University  of  Arizona 
College  of  Engineering 
Tucson,  Arizona  85721 

Professor  E.  F.  Rybicki 
The  University  of  Tulsa 
Department  of  Mechanical  Engineering 
Tulsa,  Oklahoma  74104 

Dr.  R.  Haftka 

Illinois  Institute  of  Technology 
Department  of  Mechanics  and  Mechanical 
and  Aerospace  Engineering 
Chicago,  Illinois  60616 

Professor  J.  G.  de  Oliveira 
Massachusetts  Institute  of  Technology 
Denartment  of  Ocean  Engineering 
77  Massachusetts  Avenue 
Cambridge,  Massachusetts  02139 

Dr.  Bernard  W.  Shaffer 
Polytechnic  Institute  of  New  York 
Route  110 

Farmingdale,  New  York  11735 


474:NP: 716: lab 
78u474-619 


Industry  and  Research  Institutes 

Dr.  Norman  Hobbs 
Kaaan  AviDyne 
Division  of  Kaman 
Sciences  Corporation 
Burlington,  Massachusetts  01803 

Argonne  National  Laboratory 
Library  Services  Department 
9700  South  Cass  Avenue 
Argonne,  Illinois  60440 

Dr.  M.  C.  Junger 
Cambridge  Acoustical  Associates 
54  Rindge  Avenue  Extension 
Cambridge,  Massachusetts  02140 

Mr.  J.  H.  Torrance 
General  Dynamics  Corporation 
Electric  Boat  Division 
Groton,  Connecticut  06340 

Dr.  J.  E.  Greenspon 
J.  G.  Engineering  Research  Associates 
3831  Menlo  Drive 
Baltimore,  Maryland  21215 

Newport  News  Shipbuilding  and 
Dry  Dock  Company 
Library 

Newport  News,  Virginia  23607 

Dr.  W.  F.  Bozich 

McDonnell  Douglas  Corporation 

5301  Bolsa  Avenue 

Huntington  Beach,  California  92647 

Dr.  H.  N.  Abramson 
Southwest  Research  Institute 
8500  Culebra  Road 
San  Antonio,  Texas  78284 

Dr.  R.  C.  DeHart 
Southwest  Research  Institute 
8500  Culebra  Road 
San  Antonio,  Texas  78284 

Dr.  M.  L.  Baron 
Weidlinger  Associates 
110  East  59th  Street 
New  York,  New  York  10022 


Industry  and  Research  Institutes  (Con't) 
Dr.  T.  L.  Geers 

Lockheed  Missiles  and  Space  Company 

3251  Hanover  Street 

Palo  Alto,  California  94304 

Mr.  William  Caywood 
Applied  Physics  Laboratory 
Johns  Hopkins  Road 
Laurel,  Maryland  20810 

Dr.  Robert  E.  Dunham 
Pacifica  Technology 
P.0.  Box  148 

Del  Mar,  California  92014 

Dr.  M.  F.  Kanninen 

Battel le  Columbus  Laboratories 

505  King  Avenue 

Columbus,  Ohio  43201 

Dr.  A.  A.  Hochrein 
Daedalean  Associates,  Inc. 

Springlake  Research  Road 
15110  Frederick  Road 
Woodbine,  Maryland  21797 

Dr.  James  W.  Jones 
Swanson  Service  Corporation 
P.0.  Box  5415 

Huntington  Beach,  California  92646 

Dr.  Robert  E.  Nickel 1 
Applied  Science  and  Technology 
3344  North  Torrey  Pines  Court 
Suite  220 

La  Jolla,  California  92037 

Dr.  Kevin  Thomas 
Westinghouse  Electric  Corp. 

Advanced  Reactors  Division 
P.  0.  Box  158 

Madison,  Pennsylvania  15663 

Dr.  H.  D.  Hibbitt 
Hibbitt  &  Karlsson,  Inc. 

132  George  M.  Cohan  Boulevard 
Providence,  Rhode  Island  02903 

Dr.  R.  D.  Mindlin 
89  Deer  Hill  Drive 
Ridgefield,  Connecticut  06877 
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Industry  and  Research  Institutes  (Con't) 

Dr.  Richard  E.  Dane 
Mega  Engineering 
11961  Tech  Road 

Silver  Spring,  Maryland  20904 

Mr.  G.  M.  Stanley 
Lockheed  Palo  Alto  Research 
Laboratory 
3251  Hanover  Street 
Palo  Alto,  California  94304 

Mr.  R.  L.  Cloud 

Robert  L.  Cloud  Associates,  Inc. 

2972  Adeline  Street 
Berkeley,  California  94703 


